Environmental Epidemiology Study Designs:
Disease Mapping
Demonstration
30 April 2025
Dr Darren Mayne
Public Health Epidemiologist
Illawarra Shoalhaven Local Health
District
Adjunct Senior Lecturer
Faculty of Medicine, School
of Public Health
Environmental
Epidemiology Study Designs — Disease Mapping Demonstration
(PUBH5125_DXMAP_DEMO) by
Darren
Mayne is licensed under
Creative
Commons Attribution-NonCommercial-ShareAlike 4.0 International
This session will demonstrate fitting the Besag, York, and Mollié (1991) conditional autoregressive (CAR) model for disease mapping and ecological regression applications involving areal data.
We will use a classic data set first reported by Clayton and Kaldor (1987) that describes the incidence of male lip cancer within Scottish administration districts from 1975 to 1980.
The data set was subsequently expanded to include the number of male population-years at-risk (Cressie, 1993) and the expected number of lip cancer cases and percentage of males employed in agriculture, fishing and forestry (AFF) (Stern & Cressie, 1999) in each district.
The map of Scottish administrative districts is a post-processed version of the Scotland.map S-Plus boundary file that was included in WinBUGS 1.4.3, a state-of-the-art 32-bit program for Bayesian inference when released on 6th August 2007 — now, not so much…
All data are available in a single ZIP archive available on the GeoDa Centre GuitHub repository and are described in Table 1 (columns 1 nd 2), along with the subset of variables used for this demonstration (columns 3 and 4)
Table 1: Data summary for Scottish lip cancer data set available from GeoDa Centre GitHub repository
| GeoDa variable | GeoDa description | PUBH5125 variable | PUBH5125 description |
|---|---|---|---|
| CODENO | Code converted to numeric (drop w prefix) | CODENO | Code converted to numeric (drop w prefix) |
| AREA | District polygon area | AREA | District polygon area |
| PERIMETER | District polygon perimeter | PERIMETER | District polygon perimeter |
| RECORD_ID | Unique ID | RECORD_ID | Unique ID |
| DISTRICT | District number 1-56 | dist_code | District number 1-56 |
| NAME | Name of districts from Cressie (1993) | dist_name | Name of districts from Cressie (1993) |
| CODE | District code from WinBugs | CODE | District code from WinBugs |
| CANCER | Lip cancer cases from Cressie (1993) | obs_ca | Lip cancer cases from Cressie (1993) |
| POP | Population years at risk from Cressie (1993) | dist_popn | Population years at risk from Cressie (1993) |
| CEXP | Expected cases from Lawson et al. (1999) | exp_ca | Expected cases from Stern & Cressie (1999) |
| AFF | Outdoor industry from Lawson et al. (1999) | pct_aff | Outdoor industry from Stern & Cressie (1999) |
These are the packages that need to be installed and attached to complete the analysis.
{ # Check for and install pacman package if missing
if(base::require(pacman, quietly = TRUE) == FALSE){install.packages("pacman")}
# list required packages
reqPackages <- c("tidyverse", # An opinionated collection of R packages designed for data science
"CARBayes", # Spatial Generalised Linear Mixed Models for Areal Unit Data
"coda", # A package for Output Analysis and Diagnostics for MCMC
"epitools", # A package for epidemiological analysis
"ggmcmc", # MCMC diagnostics with ggplot
"ggpubr", # A package for producing publication-ready plots
"ggrepel", # A package for position non-overlapping text labels with 'ggplot2'
"GGally", # A package to do pairplots
"knitr", # A general-purpose tool for dynamic report generation in R
"nimble", # A package for performing MCMC in R
"patchwork", # A package to combine plots
"rmarkdown", # A package for creating dynamic documents for R
"RColorBrewer", # A package of colour palettes
"scales", # Scale functions for visualization
"sf", # Simple features for R
"spdep") # A package to calculate neighbors)
# Load (and install if missing) other required packages
pacman::p_load(char = reqPackages, install = TRUE)
# Check required packages are loaded
pacman::p_loaded(char = reqPackages)
}This output shows the R environment of the analysis; it
includes dependencies loaded / attached by packages above.
## R version 4.3.3 (2024-02-29 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=English_Australia.utf8 LC_CTYPE=English_Australia.utf8 LC_MONETARY=English_Australia.utf8 LC_NUMERIC=C
## [5] LC_TIME=English_Australia.utf8
##
## time zone: Australia/Sydney
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] spdep_1.3-11 spData_2.3.4 sf_1.0-16 scales_1.4.0 RColorBrewer_1.1-3 rmarkdown_2.29 patchwork_1.3.0 nimble_1.3.0
## [9] knitr_1.50 GGally_2.2.1 ggrepel_0.9.5 ggpubr_0.6.0 ggmcmc_1.5.1.1 epitools_0.5-10.1 coda_0.19-4.1 CARBayes_6.1.1
## [17] Rcpp_1.0.12 MASS_7.3-60.0.1 lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4 purrr_1.0.2 readr_2.1.5
## [25] tidyr_1.3.1 tibble_3.2.1 ggplot2_3.5.1 tidyverse_2.0.0 pacman_0.5.1
##
## loaded via a namespace (and not attached):
## [1] DBI_1.2.2 deldir_2.0-4 s2_1.1.6 rlang_1.1.3 magrittr_2.0.3 e1071_1.7-14 compiler_4.3.3 png_0.1-8
## [9] vctrs_0.6.5 quantreg_5.97 pkgconfig_2.0.3 wk_0.9.1 shape_1.4.6.1 fastmap_1.1.1 backports_1.4.1 mcmc_0.9-8
## [17] leafem_0.2.3 utf8_1.2.4 pracma_2.4.4 tzdb_0.4.0 MatrixModels_0.5-3 xfun_0.52 satellite_1.0.5 glmnet_4.1-8
## [25] cachem_1.0.8 jsonlite_1.8.8 CARBayesdata_3.0 mapview_2.11.2 terra_1.7-83 broom_1.0.5 parallel_4.3.3 R6_2.5.1
## [33] bslib_0.7.0 stringi_1.8.3 car_3.1-2 boot_1.3-29 numDeriv_2016.8-1.1 jquerylib_0.1.4 iterators_1.0.14 base64enc_0.1-3
## [41] Matrix_1.6-5 splines_4.3.3 igraph_2.0.3 timechange_0.3.0 tidyselect_1.2.1 abind_1.4-5 rstudioapi_0.16.0 yaml_2.3.10
## [49] codetools_0.2-19 lattice_0.22-5 plyr_1.8.9 withr_3.0.0 evaluate_0.23 survival_3.5-8 ggstats_0.6.0 units_0.8-5
## [57] proxy_0.4-27 pillar_1.9.0 carData_3.0-5 KernSmooth_2.23-22 foreach_1.5.2 stats4_4.3.3 generics_0.1.3 sp_2.1-3
## [65] truncnorm_1.0-9 hms_1.1.3 class_7.3-22 glue_1.7.0 tools_4.3.3 SparseM_1.81 ggsignif_0.6.4 dotCall64_1.2
## [73] grid_4.3.3 MCMCpack_1.7-1 crosstalk_1.2.1 colorspace_2.1-0 raster_3.6-30 cli_3.6.2 spam_2.11-1 fansi_1.0.6
## [81] gtable_0.3.5 rstatix_0.7.2 sass_0.4.9 digest_0.6.35 classInt_0.4-10 htmlwidgets_1.6.4 farver_2.1.1 htmltools_0.5.8.1
## [89] lifecycle_1.0.4 leaflet_2.2.2In this section we extract, transform and load (ETL) the lip cancer
data. This can be done manually, but the folded code chunk show how it
can be done entirely withi R.
The following output section shows all the files contained in the
GeoDa Center ZIP file we downloaded using the folded
Code chunk. We will use the
scotlip/scotlip.gpkg GeoPackage data file with record index
79 at row 40 for our analysis.
# Download Scottish lip cancer data --------------------------------------------
# Available from the the GeoDa Center as a zipped archive from:
# https://geodacenter.github.io/data-and-lab/scotlip/ but can be downloaded
# directly within R
# Download ZIP file to temp directory
data.url <- "https://geodacenter.github.io/data-and-lab/data/scotlip.zip"
utils::download.file(data.url,
base::file.path(base::tempdir(), "scotlip.zip"),
mode = "wb")
# View the contents of the zip file
utils::unzip(zipfile = base::file.path(base::tempdir(), "scotlip.zip"), list = TRUE) %>%
base::as.data.frame() %>%
arrange(Name) %>%
dplyr::pull(Name)## [1] "__MACOSX/" "__MACOSX/._scotlip"
## [3] "__MACOSX/scotlip/" "__MACOSX/scotlip/._.DS_Store"
## [5] "__MACOSX/scotlip/._scotdistricts.dbf" "__MACOSX/scotlip/._scotdistricts.shp"
## [7] "__MACOSX/scotlip/._scotdistricts.shx" "__MACOSX/scotlip/._scotdistricts.txt"
## [9] "__MACOSX/scotlip/._scotlip.dbf" "__MACOSX/scotlip/._scotlip.shp"
## [11] "__MACOSX/scotlip/._scotlip.shx" "__MACOSX/scotlip/._scotlip_metadata.html"
## [13] "__MACOSX/scotlip/._scotlipdata.dbf" "__MACOSX/scotlip/._style.css"
## [15] "scotlip/" "scotlip/.DS_Store"
## [17] "scotlip/scotdistricts.csv" "scotlip/scotdistricts.dbf"
## [19] "scotlip/scotdistricts.geojson" "scotlip/scotdistricts.gpkg"
## [21] "scotlip/scotdistricts.mid" "scotlip/scotdistricts.mif"
## [23] "scotlip/scotdistricts.shp" "scotlip/scotdistricts.shx"
## [25] "scotlip/scotdistricts.sqlite" "scotlip/scotdistricts.txt"
## [27] "scotlip/scotdistricts.xlsx" "scotlip/scotlip.csv"
## [29] "scotlip/scotlip.dbf" "scotlip/scotlip.gdb/"
## [31] "scotlip/scotlip.gdb/a00000001.TablesByName.atx" "scotlip/scotlip.gdb/a00000001.gdbindexes"
## [33] "scotlip/scotlip.gdb/a00000001.gdbtable" "scotlip/scotlip.gdb/a00000001.gdbtablx"
## [35] "scotlip/scotlip.gdb/a00000002.gdbtable" "scotlip/scotlip.gdb/a00000002.gdbtablx"
## [37] "scotlip/scotlip.gdb/a00000003.gdbindexes" "scotlip/scotlip.gdb/a00000003.gdbtable"
## [39] "scotlip/scotlip.gdb/a00000003.gdbtablx" "scotlip/scotlip.gdb/a00000004.CatItemsByPhysicalName.atx"
## [41] "scotlip/scotlip.gdb/a00000004.CatItemsByType.atx" "scotlip/scotlip.gdb/a00000004.FDO_UUID.atx"
## [43] "scotlip/scotlip.gdb/a00000004.gdbindexes" "scotlip/scotlip.gdb/a00000004.gdbtable"
## [45] "scotlip/scotlip.gdb/a00000004.gdbtablx" "scotlip/scotlip.gdb/a00000004.spx"
## [47] "scotlip/scotlip.gdb/a00000005.CatItemTypesByName.atx" "scotlip/scotlip.gdb/a00000005.CatItemTypesByParentTypeID.atx"
## [49] "scotlip/scotlip.gdb/a00000005.CatItemTypesByUUID.atx" "scotlip/scotlip.gdb/a00000005.gdbindexes"
## [51] "scotlip/scotlip.gdb/a00000005.gdbtable" "scotlip/scotlip.gdb/a00000005.gdbtablx"
## [53] "scotlip/scotlip.gdb/a00000006.CatRelsByDestinationID.atx" "scotlip/scotlip.gdb/a00000006.CatRelsByOriginID.atx"
## [55] "scotlip/scotlip.gdb/a00000006.CatRelsByType.atx" "scotlip/scotlip.gdb/a00000006.FDO_UUID.atx"
## [57] "scotlip/scotlip.gdb/a00000006.gdbindexes" "scotlip/scotlip.gdb/a00000006.gdbtable"
## [59] "scotlip/scotlip.gdb/a00000006.gdbtablx" "scotlip/scotlip.gdb/a00000007.CatRelTypesByBackwardLabel.atx"
## [61] "scotlip/scotlip.gdb/a00000007.CatRelTypesByDestItemTypeID.atx" "scotlip/scotlip.gdb/a00000007.CatRelTypesByForwardLabel.atx"
## [63] "scotlip/scotlip.gdb/a00000007.CatRelTypesByName.atx" "scotlip/scotlip.gdb/a00000007.CatRelTypesByOriginItemTypeID.atx"
## [65] "scotlip/scotlip.gdb/a00000007.CatRelTypesByUUID.atx" "scotlip/scotlip.gdb/a00000007.gdbindexes"
## [67] "scotlip/scotlip.gdb/a00000007.gdbtable" "scotlip/scotlip.gdb/a00000007.gdbtablx"
## [69] "scotlip/scotlip.gdb/a00000009.gdbindexes" "scotlip/scotlip.gdb/a00000009.gdbtable"
## [71] "scotlip/scotlip.gdb/a00000009.gdbtablx" "scotlip/scotlip.gdb/a00000009.spx"
## [73] "scotlip/scotlip.gdb/a0000000a.gdbindexes" "scotlip/scotlip.gdb/a0000000a.gdbtable"
## [75] "scotlip/scotlip.gdb/a0000000a.gdbtablx" "scotlip/scotlip.gdb/a0000000a.spx"
## [77] "scotlip/scotlip.gdb/gdb" "scotlip/scotlip.gdb/timestamps"
## [79] "scotlip/scotlip.gpkg" "scotlip/scotlip.mid"
## [81] "scotlip/scotlip.mif" "scotlip/scotlip.shp"
## [83] "scotlip/scotlip.shx" "scotlip/scotlip.sqlite"
## [85] "scotlip/scotlip.xlsx" "scotlip/scotlip_metadata.html"
## [87] "scotlip/scotlipdata.dbf" "scotlip/scotlips.geojson"
## [89] "scotlip/style.css"A GeoPackage is an open-source, standards-based and platform independent file format for storing both spatial and attribute data in a single, portable file. In this code chunk, we extract the GeoPackage and then query it to identify the correct correct data layer to use.
## Extract the scotlip.gpkg geopackage data file (row 79) ----------------------
# See https://en.wikipedia.org/wiki/GeoPackage for a description of the
# geopackage data format.
# Extract the scotlip.gpkg geopackage from the ZIP file
utils::unzip(zipfile = base::file.path(base::tempdir(), "scotlip.zip"),
files = "scotlip/scotlip.gpkg",
exdir = base::tempdir(),
junkpaths = TRUE) # Only use the file name of the stored file path when extracting
# List GeoPack layers
base::cat("List GeoPackage layers:",
"\n\n",
base::as.character(base::data.frame(layers = sf::st_layers(base::file.path(base::tempdir(), "scotlip.gpkg"))$name)))As the GeoPackage contains a single layer (scotlip), we read this
into R, only to be disappointed that the wrong coordinate
reference system (WGS 80) has been assigned to the spatial data. The
correct projected coordinate reference system for these data are OSGB36 / British National Grid, which
we need to apply manually.
# Read scotlip.gpkg layer, which includes data and geometry.
# Note: I have set the coordinate reference system (CRS) to NA because the CRS
# supplied in the geopack (WGS 80) is incorrect.You can easily tell this
# by looking at he bounding box values. WGS 80 is a geographic coordinate
# system with decimal degree units, so both the minima nor maxima values
# are bounded by -180 to +180 degrees. The values here are all grater than
# 450,000, which indicate a projected coordinate system of easting and
# northing in linear meters
src_sf <- sf::st_read(dsn = base::file.path(base::tempdir(), "scotlip.gpkg"),
layer = "scotlip",
crs = NA_crs_) # Remove the wrong coordinate reference system (CRS) assigned in the geopackage## Reading layer `scotlip' from data source `C:\Users\Darren\AppData\Local\Temp\Rtmpqek4tt\scotlip.gpkg' using driver `GPKG'
## Simple feature collection with 56 features and 11 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 95631 ymin: 530297 xmax: 454570 ymax: 1203008
## CRS: NA# After some fossicking around, I was able to determine the correct projected coordinate
# system for these data is OSGB36 / British National Grid, which we now apply using
# the European Petroleum Survey Group (EPSG) code 27700
sf::st_crs(src_sf) <- 27700 # see https://epsg.io/27700 for CRS details
# Check that the CRS has been applied to the spatial data
print(src_sf)## Simple feature collection with 56 features and 11 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 95631 ymin: 530297 xmax: 454570 ymax: 1203008
## Projected CRS: OSGB36 / British National Grid
## First 10 features:
## CODENO AREA PERIMETER RECORD_ID DISTRICT NAME CODE CANCER POP CEXP AFF geom
## 1 6126 974002000 184951 1 1 Skye-Lochalsh w6126 9 28324 1.38 16 POLYGON ((214091.9 841215.2...
## 2 6016 1461990000 178224 2 2 Banff-Buchan w6016 39 231337 8.66 16 POLYGON ((383866 865862, 39...
## 3 6121 1753090000 179177 3 3 Caithness w6121 11 83190 3.04 10 POLYGON ((311487 968650, 32...
## 4 5601 898599000 128777 4 4 Berwickshire w5601 9 51710 2.53 24 POLYGON ((377180 672603, 38...
## 5 6125 5109870000 580792 5 5 Ross-Cromarty w6125 15 129271 4.26 10 POLYGON ((278680.1 882371.8...
## 6 6554 422639000 118433 6 6 Okney w6554 8 53199 2.40 24 POLYGON ((328040 1011192, 3...
## 7 6019 2267340000 259143 7 7 Moray w6019 26 245513 8.11 10 POLYGON ((348968 868724, 35...
## 8 6655 157575000 54859 8 8 Shetland w6655 7 62603 2.30 7 POLYGON ((450593 1203008, 4...
## 9 6123 4129180000 543208 9 9 Lochaber w6123 6 59183 1.98 7 POLYGON ((192646 807178, 21...
## 10 6017 2084800000 260668 10 10 Gordon w6017 20 165554 6.63 16 POLYGON ((397156 840640, 39...In this final preparation step, we wrangle the remaining attribute data (e.g. some district names need tidying up) into an intuitive analysis-ready format and enrich it by adding standardised incidence ratios and their exact Poisson 95% confidence intervals. We will also calculate indicies of estimate precision based on expected and observed lip cancer cases.
# Clean and enrich sf for data analysis
lip_sf <- src_sf %>%
# Rename columns to lower case
dplyr::rename_with(base::tolower) %>%
# Limit to selected (renamed) variables
dplyr::select(dist_code = district, dist_name = name, dist_popn = pop, obs_ca = cancer, exp_ca = cexp, pct_aff = aff) %>%
# Correct typographical errors in district names
dplyr::mutate(dplyr::across(dist_name, ~ dplyr::case_when(. == "EastKilbride" ~ "East Kilbride",
. == "EastLothian" ~ "East Lothian",
. == "NEFife" ~ "NE Fife",
. == "WesternIsles" ~ "Western Isles",
. == "WestLothian" ~ "West Lothian",
.default = .))) %>%
# Calculate standardised incidence ratios (SIR) with exact Poisson 95%
# confidence intervals using using pois.exact function from epitools.
dplyr::bind_cols(base::data.frame(epitools::pois.exact(pull(., obs_ca), pull(., exp_ca)) %>%
dplyr::select(obs_sir_est = rate,
obs_sir_l95 = lower,
obs_sir_u95 = upper))) %>%
# Calculate probability values for SIR
dplyr::mutate(obs_sir_sig = stats::dpois(obs_ca, exp_ca),
obs_sir_cat = dplyr::case_when(obs_sir_est < 1 & obs_sir_sig < 0.005 ~ "--",
obs_sir_est < 1 & obs_sir_sig < 0.025 ~ "-",
obs_sir_est > 1 & obs_sir_sig < 0.005 ~ "++",
obs_sir_est > 1 & obs_sir_sig < 0.025 ~ "+",
.default = "")) %>%
# Calculate standard and relative errors for SIR
dplyr::mutate(obs_sir_se = base::sqrt(obs_ca) / exp_ca,
obs_sir_rse = obs_sir_se / obs_sir_est, # Also estimated by RSE = 1 / base::sqrt(obs_ca),
exp_sir_rse = 1 / exp_ca) %>%
# # Add district centroids (note use of sf::st_geometry - avoids warnings about assumed attribute distributions
#
# dplyr::bind_cols(sf::st_coordinates(sf::st_centroid(sf::st_geometry(.))) %>%
# base::as.data.frame() %>%
# dplyr::rename(centroid_x = X, centroid_y = Y)) %>%
# Add variable labels
labelled::set_variable_labels(dist_code = "District code",
dist_name = "District name",
dist_popn = "Male population years at risk (1975-1980)",
obs_ca = "Observed male lip cancer cases (1975-1980)",
exp_ca = "Expected male lip cancer cases (1975-1980)",
pct_aff = "Percent of male population working in agriculture, fishing or forestry",
obs_sir_est = "Observed SIR estimate",
obs_sir_l95 = "Observed SIR lower 95% CI",
obs_sir_u95 = "Observed SIR upper 95% CI",
obs_sir_sig = "Probability value for observed SIR estimate",
obs_sir_cat = "Probability category for observed SIR estimate",
obs_sir_se = "Observed SIR standard error",
obs_sir_rse = "Observed SIR relative standard error",
exp_sir_rse = "Proportional precision of SIR from expected cases",
# centroid_x = "Centroid X coordinate (EPSG:27700)",
# centroid_y = "Centroid Y coordinate (EPSG:27700)",
geom = "Geometry (OSGB 1936 / British National Grid)") %>%
# Move geometry column to end of data set
dplyr::select(dplyr::everything(), -geom, geom)
# View the data set
print(lip_sf, digits = 2, n = 56)## Simple feature collection with 56 features and 14 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 95631 ymin: 530297 xmax: 454570 ymax: 1203008
## Projected CRS: OSGB36 / British National Grid
## dist_code dist_name dist_popn obs_ca exp_ca pct_aff obs_sir_est obs_sir_l95 obs_sir_u95 obs_sir_sig obs_sir_cat obs_sir_se obs_sir_rse exp_sir_rse
## 1 1 Skye-Lochalsh 28324 9 1.4 16 6.52 2.9822 12.38 1.3e-05 ++ 2.174 0.33 0.725
## 2 2 Banff-Buchan 231337 39 8.7 16 4.50 3.2024 6.16 3.1e-14 ++ 0.721 0.16 0.115
## 3 3 Caithness 83190 11 3.0 10 3.62 1.8063 6.47 2.5e-04 ++ 1.091 0.30 0.329
## 4 4 Berwickshire 51710 9 2.5 24 3.56 1.6266 6.75 9.3e-04 ++ 1.186 0.33 0.395
## 5 5 Ross-Cromarty 129271 15 4.3 10 3.52 1.9707 5.81 3.0e-05 ++ 0.909 0.26 0.235
## 6 6 Okney 53199 8 2.4 24 3.33 1.4391 6.57 2.5e-03 ++ 1.179 0.35 0.417
## 7 7 Moray 245513 26 8.1 10 3.21 2.0942 4.70 3.2e-07 ++ 0.629 0.20 0.123
## 8 8 Shetland 62603 7 2.3 7 3.04 1.2236 6.27 6.8e-03 + 1.150 0.38 0.435
## 9 9 Lochaber 59183 6 2.0 7 3.03 1.1121 6.60 1.2e-02 + 1.237 0.41 0.505
## 10 10 Gordon 165554 20 6.6 16 3.02 1.8426 4.66 1.5e-05 ++ 0.675 0.22 0.151
## 11 11 Western Isles 87815 13 4.4 7 2.95 1.5732 5.05 4.6e-04 ++ 0.819 0.28 0.227
## 12 12 Sutherland 37521 5 1.8 16 2.79 0.9070 6.52 2.6e-02 1.249 0.45 0.559
## 13 13 Nairn 29374 3 1.1 10 2.78 0.5728 8.12 7.1e-02 1.604 0.58 0.926
## 14 14 Wigtown 86444 8 3.3 24 2.42 1.0435 4.76 1.3e-02 + 0.855 0.35 0.302
## 15 15 NE Fife 185472 17 7.8 7 2.17 1.2632 3.47 1.8e-03 ++ 0.526 0.24 0.128
## 16 16 Kincardine 111665 9 4.5 16 1.98 0.9045 3.75 2.4e-02 + 0.659 0.33 0.220
## 17 17 Badenoch 27075 2 1.1 10 1.87 0.2264 6.75 2.0e-01 1.322 0.71 0.935
## 18 18 Ettrick 94145 7 4.2 7 1.67 0.6733 3.45 6.8e-02 0.633 0.38 0.239
## 19 19 Inverness 162867 9 5.5 7 1.63 0.7442 3.09 5.3e-02 0.542 0.33 0.181
## 20 20 Roxburgh 102697 7 4.4 10 1.58 0.6339 3.25 8.0e-02 0.596 0.38 0.225
## 21 21 Angus 263205 16 10.5 7 1.53 0.8743 2.48 2.8e-02 0.382 0.25 0.096
## 22 22 Aberdeen 583327 31 22.7 16 1.37 0.9291 1.94 1.8e-02 + 0.246 0.18 0.044
## 23 23 Argyll-Bute 190816 11 8.8 10 1.25 0.6261 2.24 9.2e-02 0.378 0.30 0.114
## 24 24 Clydesdale 163818 7 5.6 7 1.25 0.5008 2.57 1.3e-01 0.471 0.38 0.178
## 25 25 Kirkcaldy 432132 19 15.5 1 1.23 0.7394 1.92 6.3e-02 0.282 0.23 0.065
## 26 26 Dunfermline 378946 15 12.5 1 1.20 0.6722 1.98 8.1e-02 0.310 0.26 0.080
## 27 27 Nithsdale 163703 7 6.0 7 1.16 0.4660 2.39 1.4e-01 0.438 0.38 0.166
## 28 28 East Lothian 231185 10 9.0 7 1.12 0.5352 2.05 1.2e-01 0.353 0.32 0.112
## 29 29 Perth-Kinross 346041 16 14.4 10 1.11 0.6364 1.81 9.1e-02 0.278 0.25 0.070
## 30 30 West Lothian 382702 11 10.2 10 1.08 0.5383 1.93 1.2e-01 0.325 0.30 0.098
## 31 31 Cumnock-Doon 139148 5 4.8 7 1.05 0.3418 2.46 1.7e-01 0.471 0.45 0.211
## 32 32 Stewartry 65448 3 2.9 24 1.04 0.2148 3.04 2.2e-01 0.601 0.58 0.347
## 33 33 Midlothian 249667 7 7.0 10 1.00 0.4003 2.05 1.5e-01 0.376 0.38 0.142
## 34 34 Stirling 233125 8 8.5 7 0.94 0.4049 1.85 1.4e-01 0.332 0.35 0.117
## 35 35 Kyle-Carrick 319316 11 12.3 7 0.89 0.4457 1.60 1.1e-01 0.269 0.30 0.081
## 36 36 Inverclyde 296238 9 10.1 0 0.89 0.4075 1.69 1.2e-01 0.297 0.33 0.099
## 37 37 Cunninghame 391513 11 12.7 10 0.87 0.4331 1.55 1.1e-01 0.262 0.30 0.079
## 38 38 Monklands 319072 8 9.3 1 0.86 0.3694 1.69 1.3e-01 0.303 0.35 0.107
## 39 39 Dumbarton 231227 6 7.2 16 0.83 0.3058 1.81 1.4e-01 0.340 0.41 0.139
## 40 40 Clydebank 156924 4 5.3 0 0.76 0.2068 1.94 1.7e-01 0.380 0.50 0.190
## 41 41 Renfrew 617413 10 18.8 1 0.53 0.2556 0.98 1.1e-02 - 0.169 0.32 0.053
## 42 42 Falkirk 426519 8 15.8 16 0.51 0.2189 1.00 1.3e-02 - 0.179 0.35 0.063
## 43 43 Clackmannan 141294 2 4.3 16 0.46 0.0561 1.67 1.2e-01 0.327 0.71 0.231
## 44 44 Motherwell 449231 6 14.6 0 0.41 0.1505 0.89 6.0e-03 - 0.167 0.41 0.068
## 45 45 Edinburgh 1287561 19 50.7 1 0.37 0.2255 0.58 1.9e-07 -- 0.086 0.23 0.020
## 46 46 Kilmarnock 238170 3 8.2 7 0.37 0.0754 1.07 2.5e-02 0.211 0.58 0.122
## 47 47 East Kilbride 246849 2 5.6 1 0.36 0.0433 1.29 5.8e-02 0.253 0.71 0.179
## 48 48 Hamilton 312103 3 9.3 1 0.32 0.0662 0.94 1.2e-02 - 0.185 0.58 0.107
## 49 49 Glasgow 2316353 28 88.7 0 0.32 0.2099 0.46 3.5e-14 -- 0.060 0.19 0.011
## 50 50 Dundee 547016 6 19.6 1 0.31 0.1122 0.67 2.4e-04 -- 0.125 0.41 0.051
## 51 51 Cumbernauld 179194 1 3.4 1 0.29 0.0074 1.62 1.1e-01 0.291 1.00 0.291
## 52 52 Bearsden 110707 1 3.6 0 0.28 0.0070 1.54 9.7e-02 0.276 1.00 0.276
## 53 53 Eastwood 146112 1 5.7 1 0.17 0.0044 0.97 1.8e-02 - 0.174 1.00 0.174
## 54 54 Strathkelvin 246744 1 7.0 1 0.14 0.0036 0.79 6.2e-03 - 0.142 1.00 0.142
## 55 55 Tweeddale 38704 0 4.2 16 0.00 0.0000 0.89 1.6e-02 - 0.000 NaN 0.240
## 56 56 Annandale 103412 0 1.8 10 0.00 0.0000 2.10 1.7e-01 0.000 NaN 0.568
## geom
## 1 POLYGON ((214092 841215, 21...
## 2 POLYGON ((383866 865862, 4e...
## 3 POLYGON ((311487 968650, 32...
## 4 POLYGON ((377180 672603, 38...
## 5 POLYGON ((278680 882372, 29...
## 6 POLYGON ((328040 1e+06, 321...
## 7 POLYGON ((348968 868724, 35...
## 8 POLYGON ((450593 1203008, 4...
## 9 POLYGON ((192646 807178, 21...
## 10 POLYGON ((4e+05 840640, 4e+...
## 11 POLYGON ((123439 905891, 12...
## 12 POLYGON ((3e+05 965593, 291...
## 13 POLYGON ((294866 861442, 3e...
## 14 POLYGON ((246670 587314, 25...
## 15 POLYGON ((333001 723320, 34...
## 16 POLYGON ((378244 8e+05, 392...
## 17 POLYGON ((289113 829193, 3e...
## 18 POLYGON ((347058 658756, 35...
## 19 POLYGON ((232760 848360, 23...
## 20 POLYGON ((366182 640000, 37...
## 21 POLYGON ((372557 762318, 36...
## 22 POLYGON ((393827 809112, 39...
## 23 POLYGON ((209543 742400, 23...
## 24 POLYGON ((290717 655664, 30...
## 25 POLYGON ((340000 7e+05, 328...
## 26 POLYGON ((315000 7e+05, 318...
## 27 POLYGON ((265992 6e+05, 274...
## 28 POLYGON ((332787 673159, 34...
## 29 POLYGON ((316441 776777, 31...
## 30 POLYGON ((305532 679799, 31...
## 31 POLYGON ((258712 632155, 27...
## 32 POLYGON ((248144 592429, 26...
## 33 POLYGON ((314305 659828, 33...
## 34 POLYGON ((258373 734511, 26...
## 35 POLYGON ((232714 633968, 23...
## 36 POLYGON ((236768 673889, 23...
## 37 POLYGON ((227878 665244, 24...
## 38 POLYGON ((272040 672015, 28...
## 39 POLYGON ((233418 712863, 24...
## 40 POLYGON ((250189 678029, 25...
## 41 POLYGON ((247990 667816, 25...
## 42 POLYGON ((287329 691739, 29...
## 43 POLYGON ((286574 7e+05, 3e+...
## 44 POLYGON ((287710 665963, 29...
## 45 POLYGON ((311189 678784, 32...
## 46 POLYGON ((241202 652965, 24...
## 47 POLYGON ((257879 657031, 26...
## 48 POLYGON ((267205 662028, 26...
## 49 POLYGON ((251399 672324, 25...
## 50 POLYGON ((351626 732228, 33...
## 51 POLYGON ((274209 683508, 28...
## 52 POLYGON ((257360 676248, 25...
## 53 POLYGON ((252136 658209, 25...
## 54 POLYGON ((265911 681541, 27...
## 55 POLYGON ((314305 659828, 33...
## 56 POLYGON ((316714 617492, 31...Just before we jump into the analysis, I will create three additional resources:
# Dissolve (union) over the district geometries in the lip_sf sf object to obtain
# an outline (country) boundary for Scotland
sct_sf <- lip_sf %>%
dplyr::mutate(country = "Scotland") %>%
dplyr::summarise(geom = sf::st_union(geom),
.by = country)
base::print(sct_sf)## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 95631 ymin: 530297 xmax: 454570 ymax: 1203008
## Projected CRS: OSGB36 / British National Grid
## country geom
## 1 Scotland MULTIPOLYGON (((247093 5498...# Create a vector of district names to exclude when labeling maps
exc_labels <- c("Eastwood",
"Kilmarnock",
"Monklands",
"Motherwell",
"Renfrew",
"East Kilbride",
"Cumbernauld",
"Kirkcaldy",
"West Lothian",
"Inverclyde",
"Strathkelvin",
"Bearsden",
"East Lothian",
"Midlothian")# Visualization colour palette
col_pal <- c("low" = RColorBrewer::brewer.pal(11, "RdYlBu")[11],
"mid" = RColorBrewer::brewer.pal(11, "RdYlBu")[6],
"high" = RColorBrewer::brewer.pal(11, "RdYlBu")[1])
# Create a colour ramp based on colour palette
rampColours <- grDevices::colorRampPalette(col_pal)
# Visualise the colour palette
graphics::par(mar = c(0, 0, 0, 0))
graphics::image(base::seq(1, 1000), 1, base::matrix(seq(1,1000), 1000,1),
col = rampColours(1000),
pin = c(1, 5),
axes = FALSE,
ann = FALSE)### Map observed SIR -----------------------------------------------------------
# Find the absolute maximum SIR on the log scale
max_log_sir <- base::as.integer(base::max(base::abs(base::log(lip_sf[lip_sf$obs_ca != 0, ]$obs_sir_est))) / 1) + 1
# Produce the raw SIR map ------------------------------------------------------
map_obs_sir <- ggplot2::ggplot() +
# Add districts with SIR > 0
ggplot2::geom_sf(data = lip_sf %>% dplyr::filter(obs_sir_est > 0),
mapping = ggplot2::aes(fill = obs_sir_est),
col = "grey50") +
# Add districts with SIR = 0
ggplot2::geom_sf(data = lip_sf %>% dplyr::filter(obs_sir_est == 0),
mapping = ggplot2::aes(alpha = base::factor("SIR = 0", ordered = TRUE)),
fill = "grey90",
colour = "grey50",
size = 0.1) +
# Add SIR probability category
ggplot2::geom_point(data = lip_sf %>% dplyr::filter(obs_sir_cat != ""),
mapping = ggplot2::aes(shape = obs_sir_cat,
geometry = geom),
fill = "black",
stat = "sf_coordinates") +
# Add Scotland boundary
ggplot2::geom_sf(data = sct_sf,
fill = NA,
colour = "Black",
linewidth = 1) +
# Add fill aesthetic gradient mapping
ggplot2::scale_fill_gradient2("Observed SIR",
low = col_pal["low"], # blue
mid = col_pal["mid"], # yellow
high = col_pal["high"], # darkblue
midpoint = 0,
na.value = NA,
limits = c(base::exp(-max_log_sir), base::exp(max_log_sir)),
breaks = base::exp(base::seq(-max_log_sir, max_log_sir, 0.5)),
labels = base::format(base::round(exp(seq(-max_log_sir, max_log_sir, 0.5)), 2), nsmall = 2),
trans = "log") +
ggplot2::scale_alpha_manual(values = c("SIR = 0" = 1)) +
# Add shape aesthic mapping
ggplot2::scale_shape_manual("Probability value",
values = c("++" = 24,
"+" = 2,
"-" = 6,
"--" = 25),
breaks = c("++", "+", "-", "--"),
labels = c("++" = "SIR > 1, p < 0.01",
"+" = "SIR > 1, p < 0.05",
"-" = "SIR < 1, p < 0.05",
"--" = "SIR < 1, p < 0.01")) +
# Add code labels to identify districts
# ggrepel::geom_text_repel(data = lip_sf %>% dplyr::filter(! dist_name %in% exc_labels),
# inherit.aes = FALSE,
# mapping = ggplot2::aes(label = dist_name,
# geometry = geom),
# stat = "sf_coordinates",
# force = 0,
# color = "white",
# bg.color = "grey30",
# bg.r = 0.1,
# size = 9/.pt) +
# Pretty up the legend
ggplot2::guides(fill = ggplot2::guide_colourbar(frame.colour = "black",
frame.linewidth = 0.25,
ticks.colour = "black",
ticks.linewidth = 0.25,
order = 1),
alpha = ggplot2::guide_legend(title = NULL, order = 2)) +
# Pretty up the map
ggplot2::theme_bw() +
ggplot2::theme(axis.title = ggplot2::element_blank(),
axis.text = ggplot2::element_blank(),
axis.ticks = ggplot2::element_blank(),
legend.position = "inside",
legend.margin = ggplot2::margin(-3, 0, 0, 0),
legend.position.inside = c(0.15, 0.85),
panel.grid = ggplot2::element_blank(),
text = ggplot2::element_text(size = 12))
# Map relative standard errors for the observed SIR ----------------------------
map_obs_rse <- ggplot2::ggplot(data = lip_sf) +
# Add relative standard errors for districts with SIR > 0
# Add districts with SIR > 0
ggplot2::geom_sf(data = lip_sf %>% dplyr::filter(obs_sir_est > 0),
mapping = ggplot2::aes(fill = obs_sir_rse),
col = "grey50") +
# Add districts with SIR = 0
ggplot2::geom_sf(data = lip_sf %>% dplyr::filter(obs_sir_est == 0),
mapping = ggplot2::aes(alpha = base::factor("SIR = 0", ordered = TRUE)),
fill = "grey90",
colour = "grey50",
size = 0.1) +
# Add Scotland boundary
ggplot2::geom_sf(data = sct_sf,
fill = NA,
colour = "Black",
linewidth = 1) +
# Add code labels to identify districts
ggrepel::geom_text_repel(data = (lip_sf %>%
dplyr::filter(dist_name %in% c("Berwickshire",
"Inverness",
"Aberdeen",
"Glasgow"))),
inherit.aes = FALSE,
mapping = ggplot2::aes(label = dist_name,
geometry = geom),
stat = "sf_coordinates",
force = 0,
color = "white",
bg.color = "grey30",
bg.r = 0.1,
size = 12/.pt) +
# Add fill aesthetic gradient mapping
ggplot2::scale_fill_gradient2("Relative standard error",
low = col_pal["low"], # blue
mid = col_pal["mid"], # yellow
high = col_pal["high"], # darkblue
midpoint = 0.5,
na.value = NA,
limits = c(0, 1),
breaks = base::seq(0, 1, 0.1),
labels = base::seq(0, 1, 0.1)) +
ggplot2::scale_alpha_manual(values = c("SIR = 0" = 1)) +
ggplot2::guides(fill = ggplot2::guide_colourbar(frame.colour = "black",
frame.linewidth = 0.25,
ticks.colour = "black",
ticks.linewidth = 0.25,
order = 1),
alpha = ggplot2::guide_legend(title = NULL, order = 2)) +
ggplot2::theme_bw() +
ggplot2::theme(axis.title = ggplot2::element_blank(),
axis.text = ggplot2::element_blank(),
axis.ticks = ggplot2::element_blank(),
legend.position = "inside",
legend.margin = ggplot2::margin(-3, 0, 0, 0),
# legend.position.inside = c(0.15, 0.88),
legend.position.inside = c(0.18, 0.90),
legend.text = ggplot2::element_text(size = 10),
panel.grid = ggplot2::element_blank(),
text = ggplot2::element_text(size = 12))
# Map expected precision of SIR proportional to expected cancer cases -----------
map_exp_rse <- ggplot2::ggplot(data = lip_sf) +
# Add relative standard errors for districts with SIR > 0
# Add districts with E(R) > 0
ggplot2::geom_sf(data = lip_sf, # %>% dplyr::filter(obs_sir_est > 0),
mapping = ggplot2::aes(fill = exp_sir_rse),
col = "grey50") +
# Add Scotland boundary
ggplot2::geom_sf(data = sct_sf,
fill = NA,
colour = "Black",
linewidth = 1) +
# Add fill aesthetic gradient mapping
ggplot2::scale_fill_gradient2(base::expression("SE(R)" %prop% "1/E"), # "Expected precision of SIR",
low = col_pal["low"], # blue
mid = col_pal["mid"], # yellow
high = col_pal["high"], # darkblue
midpoint = 0.5,
na.value = NA,
limits = c(0, 1),
breaks = base::seq(0, 1, 0.1),
labels = base::seq(0, 1, 0.1)) +
# ggplot2::scale_alpha_manual(values = c("SIR = 0" = 1)) +
ggplot2::guides(fill = ggplot2::guide_colourbar(frame.colour = "black",
frame.linewidth = 0.25,
ticks.colour = "black",
ticks.linewidth = 0.25,
order = 1),
alpha = ggplot2::guide_legend(title = NULL, order = 2)) +
ggplot2::theme_bw() +
ggplot2::theme(axis.title = ggplot2::element_blank(),
axis.text = ggplot2::element_blank(),
axis.ticks = ggplot2::element_blank(),
legend.position = "inside",
legend.margin = ggplot2::margin(-3, 0, 0, 0),
# legend.position.inside = c(0.15, 0.88),
legend.position.inside = c(0.13, 0.905),
legend.text = ggplot2::element_text(size = 10),
panel.grid = ggplot2::element_blank(),
text = ggplot2::element_text(size = 12))
## Output observed and BYM smoothed SIR for lip cancer risk --------------------
(map_obs_sir | map_obs_rse | map_exp_rse) +
patchwork::plot_annotation(title = "Observed SIR estimates and p-values and observed and expected relative standard errors for male lip cancer in Scotland (1975-1980)",
theme = theme(plot.title = element_text(hjust = 0.5)))Prior to conducting the BYM CAR analysis, we need to construct the
queen adjacency matrix. This is easily using the
spdep::poly2nb function.
### Create Queen adjacency (contiguity) list for each district -----------------
# The spdep::poly2nb function constructs a list of neighbours from a polygon object
adj_nb_list <- spdep::poly2nb(lip_sf,
queen = TRUE,
row.names = base::row.names(lip_sf))## Warning in spdep::poly2nb(lip_sf, queen = TRUE, row.names = base::row.names(lip_sf)): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.
## Warning in spdep::poly2nb(lip_sf, queen = TRUE, row.names = base::row.names(lip_sf)): neighbour object has 4 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.
## Neighbour list object:
## Number of regions: 56
## Number of nonzero links: 234
## Percentage nonzero weights: 7.461735
## Average number of links: 4.178571
## 3 regions with no links:
## 6, 8, 11
## 4 disjoint connected subgraphs## [[1]]
## [1] 5 9 19
##
## [[2]]
## [1] 7 10
##
## [[3]]
## [1] 12
##
## [[4]]
## [1] 18 20 28
##
## [[5]]
## [1] 1 12 19
##
## [[6]]
## [1] 0The warning message is telling us three districts (districts 6, 8 and 11) have no neighbours, which is reasonable as there are three island areas (districts), and that one has (district 4) has disjoint connected subgraphs. This likely means the area is sharing multiple borders (edges) with the same neighbour. This can occur when there a “slivers” — small gaps — in your spatial data between bordering spatial units. Given the map available for Scotland looks like a five-year old has cut it out of coloured paper with blunt craft scissors and from memory, this almost certainly likely. We can check these warnings are OK via some additional maps.
# Provide district map for reference -------------------------------------------
adj_dist_map <- ggplot2::ggplot() +
# Add districts with SIR > 0
ggplot2::geom_sf(data = lip_sf,
fill = NA,
col = "black",
linewidth = 0.25) +
# Add Scotland boundary
ggplot2::geom_sf(data = sct_sf,
fill = NA,
colour = "Black",
linewidth = 1) +
# Add code labels to identify districts
ggplot2::geom_text(data = lip_sf,
inherit.aes = FALSE,
mapping = ggplot2::aes(label = dist_code,
geometry = geom),
stat = "sf_coordinates",
color = "black",
size = 12/.pt) +
ggplot2::theme_bw() +
ggplot2::theme(axis.title = ggplot2::element_blank(),
axis.text = ggplot2::element_blank(),
axis.ticks = ggplot2::element_blank(),
legend.position = "inside",
legend.margin = ggplot2::margin(-3, 0, 0, 0),
# legend.position.inside = c(0.15, 0.88),
legend.position.inside = c(0.13, 0.905),
legend.text = ggplot2::element_text(size = 10),
panel.grid = ggplot2::element_blank(),
text = ggplot2::element_text(size = 12))
# Map some neighbours ----------------------------------------------------------
adj_nb_map <- ggplot2::ggplot() +
# Add districts
ggplot2::geom_sf(data = lip_sf,
colour = "grey50",
fill = NA,
linewidth = 0.1) +
# Areas of interest
# 4 = Berwickshire
# 19 = Inverness
# 22 = Aberdeen
# 49 = Glasgow
ggplot2::geom_sf(data = lip_sf %>% filter(row.names(.) %in% c(4, 19, 22, 49)),
mapping = ggplot2::aes(fill = "Analysis district"),
colour = "black",
linewidth = 0.1) +
# Neighbours of areas of interest
ggplot2::geom_sf(data = lip_sf %>% filter(row.names(.) %in% unlist(adj_nb_list[c(4, 19, 22, 49)])),
mapping = aes(fill = "Queen neighbours"),
colour = "black",
linewidth = 0.1) +
# Add Scotland boundary
ggplot2::geom_sf(data = sct_sf,
fill = NA,
colour = "Black",
linewidth = 1) +
# Add code labels to identify districts
ggrepel::geom_text_repel(data = (lip_sf %>%
dplyr::filter(dist_name %in% c("Berwickshire",
"Inverness",
"Aberdeen",
"Glasgow"))),
inherit.aes = FALSE,
mapping = ggplot2::aes(label = dist_name,
geometry = geom),
stat = "sf_coordinates",
force = 0,
color = "white",
bg.color = "grey30",
bg.r = 0.1,
size = 12/.pt) +
# Add fill aesthetic gradient mapping
ggplot2::scale_fill_manual("Adjacency examples",
breaks = c("Analysis district",
"Queen neighbours"),
values = c(rgb(230, 70, 38, maxColorValue = 255),
rgb(255, 212, 102, maxColorValue = 255))) +
# Add X and Y coordinate labels
ggplot2::theme_bw() +
ggplot2::theme(axis.title = ggplot2::element_blank(),
axis.text = ggplot2::element_blank(),
axis.ticks = ggplot2::element_blank(),
legend.position = "inside",
legend.margin = ggplot2::margin(0, 0, 0, 0),
legend.position.inside = c(0.2, 0.93),
panel.grid = ggplot2::element_blank(),
text = ggplot2::element_text(size = 12))
# Plot neighbours connectivity graph -------------------------------------------
# Genereate connectivity arcs between spatial units and their neighbours
# From: https://mbjoseph.github.io/posts/2018-12-27-plotting-spatial-neighbors-in-ggplot2/
# The spdep::nb2lines function creates linestring geometries linking the centroids
# of neighbouring districts
adj_nb_arc <- spdep::nb2lines(adj_nb_list,
coords = st_centroid(st_geometry(lip_sf)))
# Calculate the number of neighbours for each district by summarising the length
# of each element in the adj_nb_list
adj_nb_num <- lapply(adj_nb_list, length) %>%
do.call(rbind, .) %>%
as.data.frame() %>%
setNames("neighbours")
## Map adjacency marks and symbolise by number of neighbours -------------------
adj_con_map <- ggplot() +
# Add districts
geom_sf(data = lip_sf %>% bind_cols(adj_nb_num),
aes(fill = neighbours),
colour = "grey50",
linewidth = 0.1,
alpha = 0.75) +
# Add neighbour connectivity arcs
geom_sf(data = adj_nb_arc,
colour = "black",
linewidth = 0.1) +
# Add district centroids
geom_sf(data = st_centroid(st_geometry(lip_sf)),
shape = 16,
colour = "black",
size = 2.75) +
# Add code labels to identify districts
ggrepel::geom_text_repel(data = lip_sf %>% dplyr::filter(! dist_name %in% exc_labels),
inherit.aes = FALSE,
mapping = ggplot2::aes(label = dist_name,
geometry = geom),
stat = "sf_coordinates",
force = 0,
color = "white",
bg.color = "grey30",
bg.r = 0.1,
size = 12/.pt) +
# Add fill aesthetic gradient mapping
scale_fill_gradient("Neighbours",
low = "lightyellow",
high = "darkred",
limits = c(0, 11),
breaks = seq(0, 11, 1),
# labels = c("1", "", "", "", "", "6", "", "", "", "", "", "12"),
guide = guide_colorbar(frame.colour = "black",
frame.linewidth = 0.25,
ticks.colour = "black",
ticks.linewidth = 0.25)) +
theme_bw() +
theme(axis.title = ggplot2::element_blank(),
axis.text = ggplot2::element_blank(),
axis.ticks = ggplot2::element_blank(),
legend.position = "inside",
legend.margin = ggplot2::margin(0, 0, 0, 0),
legend.position.inside = c(0.15, 0.90),
panel.grid = ggplot2::element_blank(),
text = ggplot2::element_text(size = 12))
## Output adjacency maps -------------------------------------------------------
(adj_dist_map | adj_nb_map | adj_con_map) +
patchwork::plot_annotation(title = "Spatial neighbours and links maps for Scottish administrative districts",
theme = theme(plot.title = element_text(hjust = 0.5)))Based on these maps, we can be confident the three districts without
neighbours are island areas (Western Isles, Orkney and Shetland) and
that district 4 (Berwickshire) is connected to the correct number of
neighbours (n=3). We can now convert the spdep adjacency
matrix to WinBUGS weight matrix format, which is required for fitting
BYM CAR models in nimble.
### Output adjacency matrix in WinBUGS format for use in BYM models ------------
# The spdep::nb2WB function can be used to make WinBUGS neighbours list
adj_nb_wgt <- spdep::nb2WB(nb = adj_nb_list)
# View WinBUIGS adjacency list
base::print(adj_nb_wgt)## $adj
## [1] 5 9 19 7 10 12 18 20 28 1 12 19 2 10 13 16 17 1 17 19 23 29 2 7 16 22 3 5 7 17 19 31 32 35 25 29 7 10 17 21 22 29 7 9 13 16 19 29 4 20 28 33 55 56 1
## [56] 5 9 13 17 4 18 56 16 29 50 10 16 9 29 34 39 27 30 31 44 47 48 55 56 15 26 29 25 29 43 24 31 32 56 4 18 33 45 9 15 16 17 21 23 25 26 34 43 50 24 38 42 44 45 55
## [111] 14 24 27 32 35 46 47 14 27 31 18 28 45 55 23 29 39 40 42 43 51 52 54 14 31 37 46 37 41 35 36 41 46 30 42 44 49 51 54 23 34 40 34 39 49 52 36 37 46 49 53 30 34 38 43
## [166] 51 26 29 34 42 24 30 38 48 49 28 30 33 55 31 35 37 41 47 53 24 31 46 48 49 53 24 44 47 49 38 40 41 44 47 48 52 53 54 21 29 34 38 42 54 34 40 49 54 41 46 47 49 34 38
## [221] 49 51 52 18 24 30 33 45 56 18 20 24 27 55
##
## $weights
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [83] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [165] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##
## $num
## [1] 3 2 1 3 3 0 5 0 5 4 0 2 3 3 2 6 6 6 5 3 3 2 4 8 3 3 4 4 11 6 7 3 4 9 4 2 4 6 3 4 5 5 4 5 4 6 6 4 9 2 4 4 4 5 6
## [56] 5# The *adj* vector defines the district indexes for the neighbours of the ith
# analysis unit; the *weights* vector indicates the weight assigned to each
# neighbour when calculating the conditional autoregression spatial term
# in the BYM model; and *num* indicates how many of the neighbours in *adj*
# border the ith analysis unit. For example, the first element of *num* indicates
# that the district of Skye-Lochalsh (row 1 in *lip_sf*) has three neighbours,
# which *adj* indicates correspond to row indexes 5, 9 and 19 (Ross-Cromarty,
# Lochaber and Inverness districts) in *lip_sf*, and are assigned weights of 1,
# 1 and 1 when calculating the CAR spatial term for Skye-Lochalsh.We now have everything necessary to run a BYM CAR disease mapping
model. All that is left is to format the data in the way
nimble expects it to be, write and compile our model, and
then run and summarise its output.
### Assemble objects needed to run a BYM disease model -------------------------
## Obtain the number of districts ----------------------------------------------
N <- dim(lip_sf)[1]
## Format SIR data for NIMBLE in a list ----------------------------------------
bym_data <- base::list(obs = lip_sf$obs_ca) # Observed lip cancer cases
## List constants required by NIMBLE to fit BYM model --------------------------
bym_consts <-base::list(N = N, # Number of districts
exp = lip_sf$exp_ca, # Expected lip cancer cases
# Adjacency matrix values used to specify CAR distribution
L = length(adj_nb_wgt$weights), # Total number of neighbours
adj = adj_nb_wgt$adj, # Vector of neighbours for each ditrict
num = adj_nb_wgt$num, # Vector of neighbour counts for each district
weights = adj_nb_wgt$weights) # Weight for each neighbour (all 1 for a queen matrix
base::print(bym_consts)## $N
## [1] 56
##
## $exp
## [1] 1.38 8.66 3.04 2.53 4.26 2.40 8.11 2.30 1.98 6.63 4.40 1.79 1.08 3.31 7.84 4.55 1.07 4.18 5.53 4.44 10.46 22.67 8.77 5.62 15.47 12.49 6.04
## [28] 8.96 14.37 10.20 4.75 2.88 7.03 8.53 12.32 10.10 12.68 9.35 7.20 5.27 18.76 15.78 4.32 14.63 50.72 8.20 5.59 9.34 88.66 19.62 3.44 3.62 5.74 7.03
## [55] 4.16 1.76
## attr(,"label")
## [1] "Expected male lip cancer cases (1975-1980)"
##
## $L
## [1] 234
##
## $adj
## [1] 5 9 19 7 10 12 18 20 28 1 12 19 2 10 13 16 17 1 17 19 23 29 2 7 16 22 3 5 7 17 19 31 32 35 25 29 7 10 17 21 22 29 7 9 13 16 19 29 4 20 28 33 55 56 1
## [56] 5 9 13 17 4 18 56 16 29 50 10 16 9 29 34 39 27 30 31 44 47 48 55 56 15 26 29 25 29 43 24 31 32 56 4 18 33 45 9 15 16 17 21 23 25 26 34 43 50 24 38 42 44 45 55
## [111] 14 24 27 32 35 46 47 14 27 31 18 28 45 55 23 29 39 40 42 43 51 52 54 14 31 37 46 37 41 35 36 41 46 30 42 44 49 51 54 23 34 40 34 39 49 52 36 37 46 49 53 30 34 38 43
## [166] 51 26 29 34 42 24 30 38 48 49 28 30 33 55 31 35 37 41 47 53 24 31 46 48 49 53 24 44 47 49 38 40 41 44 47 48 52 53 54 21 29 34 38 42 54 34 40 49 54 41 46 47 49 34 38
## [221] 49 51 52 18 24 30 33 45 56 18 20 24 27 55
##
## $num
## [1] 3 2 1 3 3 0 5 0 5 4 0 2 3 3 2 6 6 6 5 3 3 2 4 8 3 3 4 4 11 6 7 3 4 9 4 2 4 6 3 4 5 5 4 5 4 6 6 4 9 2 4 4 4 5 6
## [56] 5
##
## $weights
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [83] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [165] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1## Write the model specification -----------------------------------------------
# NOTE: Don't use R namespaces as models are run by the nimble C++ compiler
bym_code <- nimble::nimbleCode(
{
for (i in 1:N){
# Poisson likelihood for observed counts
obs[i] ~ dpois(lambda[i])
log(lambda[i]) <- alpha + s[i] + u[i] + log(exp[i])
# Prior for the area-specific unstructured random effect
u[i] ~ dnorm(0, tau = tau.u)
# Smoothed area-specific SIR
smoothed.sir[i] <- exp(alpha + s[i] + u[i])
# Overall residual area-specific cancer risk above or below the Scottish average
residual.sir[i] <- exp(s[i] + u[i])
# Residual area-specific cancer risk due to unobserved and spatially structured factors (s)
residual.sir.s[i] <- exp(s[i])
# Residual area-specific lip cancer risk due to unobserved and spatially unstructured factors (u)
residual.sir.u[i] <- exp(u[i])
# Posterior probabilities for residual risk - is the SIR for the ith district higher or lower than 1
residual.sir.gt1[i] <- 1 - step(1 - residual.sir[i]) # Posterior probability that SIR > 1
residual.sir.lt1[i] <- 1 - step(residual.sir[i] - 1) # Posterior probability that SIR < 1
}
# Prior for the area-specific spatially structured random effect (ICAR)
s[1:N] ~ dcar_normal(adj[1:L],
weights[1:L],
num[1:N],
tau.s,
zero_mean = 1)
# Vague uniform prior for intercept
alpha ~ dflat()
# Overall SIR for the study region (lip cancer risk common to all areas)
overall.sir <- exp(alpha)
# Hyperprior distribution on inverse variance (precision) for the unstructured variance component (u)
tau.u ~ dgamma(1, 0.01)
# Variance of u (unstructured variance component)
sigma2.u <- 1/tau.u
# Hyperprior distribution on inverse variance (precision) for the spatially structured variance component (s)
tau.s ~ dgamma(1, 0.01)
# variance of s (spatial variance component)
sigma2.s <- 1/tau.s
}
)
base::print(bym_code)## {
## for (i in 1:N) {
## obs[i] ~ dpois(lambda[i])
## log(lambda[i]) <- alpha + s[i] + u[i] + log(exp[i])
## u[i] ~ dnorm(0, tau = tau.u)
## smoothed.sir[i] <- exp(alpha + s[i] + u[i])
## residual.sir[i] <- exp(s[i] + u[i])
## residual.sir.s[i] <- exp(s[i])
## residual.sir.u[i] <- exp(u[i])
## residual.sir.gt1[i] <- 1 - step(1 - residual.sir[i])
## residual.sir.lt1[i] <- 1 - step(residual.sir[i] - 1)
## }
## s[1:N] ~ dcar_normal(adj[1:L], weights[1:L], num[1:N], tau.s,
## zero_mean = 1)
## alpha ~ dflat()
## overall.sir <- exp(alpha)
## tau.u ~ dgamma(1, 0.01)
## sigma2.u <- 1/tau.u
## tau.s ~ dgamma(1, 0.01)
## sigma2.s <- 1/tau.s
## }## Specify initial values for nodes as a data list -----------------------------
# alpha = intercept
# u = prior for unstructured random effect
# tau.u = hyperprior for the precision of u
# s = prior for structured random effect
# tau.s = hyperprior for the precision of s
bym_inits <- list(list(alpha = 0.01, # MCMC chain 1
u = base::rep(0.01, times = N),
tau.u = 10,
s = base::rep(0.01, times = N),
tau.s = 10),
list(alpha = 0.5, # MCMC chain 2
u = base::rep(-0.01, times = N),
tau.u = 1,
s = base::rep(-0.01, times = N),
tau.s = 1))
base::print(bym_inits)## [[1]]
## [[1]]$alpha
## [1] 0.01
##
## [[1]]$u
## [1] 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01
## [34] 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01
##
## [[1]]$tau.u
## [1] 10
##
## [[1]]$s
## [1] 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01
## [34] 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01
##
## [[1]]$tau.s
## [1] 10
##
##
## [[2]]
## [[2]]$alpha
## [1] 0.5
##
## [[2]]$u
## [1] -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01
## [28] -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01
## [55] -0.01 -0.01
##
## [[2]]$tau.u
## [1] 1
##
## [[2]]$s
## [1] -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01
## [28] -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01
## [55] -0.01 -0.01
##
## [[2]]$tau.s
## [1] 1## Select parameters (nodes) to monitor ----------------------------------------
bym_params <- c("alpha",
"u",
"s",
"sigma2.s",
"sigma2.u",
"overall.sir",
"residual.sir",
"residual.sir.s",
"residual.sir.u",
"smoothed.sir",
"residual.sir.gt1",
"residual.sir.lt1")
print(bym_params)## [1] "alpha" "u" "s" "sigma2.s" "sigma2.u" "overall.sir" "residual.sir" "residual.sir.s"
## [9] "residual.sir.u" "smoothed.sir" "residual.sir.gt1" "residual.sir.lt1"## Specify the MCMC sampler settings -------------------------------------------
n_chains <- 2 # Number of MCMC chains to run
n_burnin <- 50000 # Number of iterations to use as model burn-in
n_iter <- n_burnin * 2 # Number of iterations per chain
n_thin <- 1 # Thinning interval for chain, i.e., keep every nth sample
# Burn-in samples per chain (n = 2) - used to converge model
(n_burnin / n_thin) # 5000 samples used to burn the model in# Inference samples per chain (n = 2)
((n_iter - n_burnin) / n_thin) # 5000 samples from the posterior distribution summarised for inference### Fit the model --------------------------------------------------------------
# The nimble showCompilerOutput and nimbleMCMC options are set to TRUE for
# demonstration purposes only. You can exclude these options if you prefer
nimbleOptions(showCompilerOutput = TRUE)
bym_samples <- nimble::nimbleMCMC(code = bym_code,
data = bym_data,
constants = bym_consts,
inits = bym_inits,
monitors = bym_params,
niter = n_iter,
nburnin = n_burnin,
thin = n_thin,
nchains = n_chains,
setSeed = c(9, 10), # Seed for randome number generator - important for reproducibility - need one seed per chain
progressBar = TRUE,
samplesAsCodaMCMC = TRUE,
summary = TRUE,
WAIC = TRUE)## Defining model
## Building model
## Setting data and initial values
## Running calculate on model
## [Note] Any error reports that follow may simply reflect missing values in model variables.
## Checking model sizes and dimensions
## Checking model calculations
## Compiling
## [Note] This may take a minute.
## [Note] On some systems there may be some compiler warnings that can be safely ignored.
## g++ -std=gnu++17 -I"C:/Software/R/R-43~1.3/include" -DNDEBUG -DR_NO_REMAP -I"C:\Software\R\R-4.3.3\library\nimble\include" -I"" -DEIGEN_MPL2_ONLY=1 -Wno-misleading-indentation -Wno-ignored-attributes -Wno-deprecated-declarations -std=c++11 -I"C:/Software/R/rtools43/x86_64-w64-mingw32.static.posix/include" -O2 -Wall -mfpmath=sse -msse2 -mstackrealign -c dynamicRegistrations_04_30_14_18_05.cpp -o dynamicRegistrations_04_30_14_18_05.o
## g++ -std=gnu++17 -shared -s -static-libgcc -o dynamicRegistrations_04_30_14_18_05.dll tmp.def dynamicRegistrations_04_30_14_18_05.o -LC:\Software\R\R-4.3.3\library\nimble\CppCode -lnimble_x64 -lRlapack -lRblas -LC:/Software/R/rtools43/x86_64-w64-mingw32.static.posix/lib/x64 -LC:/Software/R/rtools43/x86_64-w64-mingw32.static.posix/lib -LC:/Software/R/R-43~1.3/bin/x64 -lR
## using C++ compiler: 'G__~1.EXE (GCC) 12.2.0'
## g++ -std=gnu++17 -I"C:/Software/R/R-43~1.3/include" -DNDEBUG -DR_NO_REMAP -I"C:\Software\R\R-4.3.3\library\nimble\include" -I"" -DEIGEN_MPL2_ONLY=1 -Wno-misleading-indentation -Wno-ignored-attributes -Wno-deprecated-declarations -std=c++11 -I"C:/Software/R/rtools43/x86_64-w64-mingw32.static.posix/include" -O2 -Wall -mfpmath=sse -msse2 -mstackrealign -c P_1_code_MID_1.cpp -o P_1_code_MID_1.o
## g++ -std=gnu++17 -I"C:/Software/R/R-43~1.3/include" -DNDEBUG -DR_NO_REMAP -I"C:\Software\R\R-4.3.3\library\nimble\include" -I"" -DEIGEN_MPL2_ONLY=1 -Wno-misleading-indentation -Wno-ignored-attributes -Wno-deprecated-declarations -std=c++11 -I"C:/Software/R/rtools43/x86_64-w64-mingw32.static.posix/include" -O2 -Wall -mfpmath=sse -msse2 -mstackrealign -c P_1_code_MID_1_nfCode.cpp -o P_1_code_MID_1_nfCode.o
## g++ -std=gnu++17 -shared -s -static-libgcc -o P_1_code_MID_1_04_30_14_18_05.dll tmp.def P_1_code_MID_1.o P_1_code_MID_1_nfCode.o -LC:\Software\R\R-4.3.3\library\nimble\CppCode -lnimble_x64 -lRlapack -lRblas -LC:/Software/R/rtools43/x86_64-w64-mingw32.static.posix/lib/x64 -LC:/Software/R/rtools43/x86_64-w64-mingw32.static.posix/lib -LC:/Software/R/R-43~1.3/bin/x64 -lR
## using C++ compiler: 'G__~1.EXE (GCC) 12.2.0'
## g++ -std=gnu++17 -I"C:/Software/R/R-43~1.3/include" -DNDEBUG -DR_NO_REMAP -I"C:\Software\R\R-4.3.3\library\nimble\include" -I"" -DEIGEN_MPL2_ONLY=1 -Wno-misleading-indentation -Wno-ignored-attributes -Wno-deprecated-declarations -std=c++11 -I"C:/Software/R/rtools43/x86_64-w64-mingw32.static.posix/include" -O2 -Wall -mfpmath=sse -msse2 -mstackrealign -c P_1_MCMC.cpp -o P_1_MCMC.o
## g++ -std=gnu++17 -I"C:/Software/R/R-43~1.3/include" -DNDEBUG -DR_NO_REMAP -I"C:\Software\R\R-4.3.3\library\nimble\include" -I"" -DEIGEN_MPL2_ONLY=1 -Wno-misleading-indentation -Wno-ignored-attributes -Wno-deprecated-declarations -std=c++11 -I"C:/Software/R/rtools43/x86_64-w64-mingw32.static.posix/include" -O2 -Wall -mfpmath=sse -msse2 -mstackrealign -c P_1_code_MID_1.cpp -o P_1_code_MID_1.o
## g++ -std=gnu++17 -shared -s -static-libgcc -o P_1_MCMC_04_30_14_18_29.dll tmp.def P_1_MCMC.o P_1_code_MID_1.o -LC:\Software\R\R-4.3.3\library\nimble\CppCode -lnimble_x64 -lRlapack -lRblas -LC:/Software/R/rtools43/x86_64-w64-mingw32.static.posix/lib/x64 -LC:/Software/R/rtools43/x86_64-w64-mingw32.static.posix/lib -LC:/Software/R/R-43~1.3/bin/x64 -lR
## using C++ compiler: 'G__~1.EXE (GCC) 12.2.0'
## In file included from C:\Software\R\R-4.3.3\library\nimble\include/nimble/nimbleEigen.h:162,
## from C:\Software\R\R-4.3.3\library\nimble\include/nimble/EigenTypedefs.h:33,
## from P_1_MCMC.cpp:8:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/nimbleEigenNimArr.h: In instantiation of 'void assignVectorToNimArr(NimArrOutput&, const VectorInput&) [with NimArrOutput = NimArr<1, int>; VectorInput = std::vector<int>]':
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/nimbleEigenNimArr.h:339:56: required from 'void setWhich(NimArrOutput&, const DerivedBool&) [with NimArrOutput = NimArr<1, int>; DerivedBool = Eigen::CwiseBinaryOp<Eigen::internal::scalar_cmp_op<double, double, Eigen::internal::cmp_LT>, const Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::Array<double, -1, -1, 1, -1, -1> >, const Eigen::ArrayWrapper<const Eigen::CwiseBinaryOp<Eigen::internal::scalar_quotient_op<double, double>, const Eigen::Transpose<Eigen::Block<Eigen::Map<Eigen::Matrix<double, -1, -1>, 0, Eigen::Stride<0, 0> >, -1, -1, false> >, const Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, const Eigen::Matrix<double, -1, -1, 1, -1, -1> > > > >]'
## P_1_MCMC.cpp:1265:10: required from here
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/nimbleEigenNimArr.h:318:22: warning: comparison of integer expressions of different signedness: 'int' and 'std::vector<int>::size_type' {aka 'long long unsigned int'} [-Wsign-compare]
## 318 | if(output.size() != input.size()) {
## | ~~~~~~~~~~~~~~^~~~~~~~~~~~~~~
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/nimbleEigenNimArr.h:322:22: warning: comparison of integer expressions of different signedness: 'int' and 'std::vector<int>::size_type' {aka 'long long unsigned int'} [-Wsign-compare]
## 322 | for(int i = 0; i < input.size(); i++) output[i] = input[i];
## | ~~^~~~~~~~~~~~~~
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/nimbleEigenNimArr.h:325:22: warning: comparison of integer expressions of different signedness: 'int' and 'std::vector<int>::size_type' {aka 'long long unsigned int'} [-Wsign-compare]
## 325 | for(int i = 0; i < input.size(); i++) output.valueNoMap(i) = input[i];
## | ~~^~~~~~~~~~~~~~
## In file included from C:\Software\R\R-4.3.3\library\nimble\include/nimble/EigenTypedefs.h:28:
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'nimSmartPtr<T>::~nimSmartPtr() [with T = waicNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:116:29:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '*this.nimSmartPtr<waicNimbleList>::realPtr' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'nimSmartPtr<T>::~nimSmartPtr() [with T = waicDetailsNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:116:29:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '*this.nimSmartPtr<waicDetailsNimbleList>::realPtr' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'nimSmartPtr<T>::~nimSmartPtr() [with T = waicNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:116:29,
## inlined from 'nimSmartPtr<T>::~nimSmartPtr() [with T = waicNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:117:3:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '*this.nimSmartPtr<waicNimbleList>::realPtr' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'nimSmartPtr<T>::~nimSmartPtr() [with T = waicDetailsNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:116:29,
## inlined from 'nimSmartPtr<T>::~nimSmartPtr() [with T = waicDetailsNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:117:3:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '*this.nimSmartPtr<waicDetailsNimbleList>::realPtr' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'void nimSmartPtr<T>::setPtrFromT(T*&) [with T = waicNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:55:39,
## inlined from 'void nimSmartPtr<T>::setPtrFromVoidPtr(void*&) [with T = waicNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:64:16:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '*this.nimSmartPtr<waicNimbleList>::realPtr' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'void nimSmartPtr<T>::setPtrFromT(T*&) [with T = waicDetailsNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:55:39,
## inlined from 'void nimSmartPtr<T>::setPtrFromVoidPtr(void*&) [with T = waicDetailsNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:64:16:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '*this.nimSmartPtr<waicDetailsNimbleList>::realPtr' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'void nimSmartPtr<T>::setPtr(const nimSmartPtr<T>&) [with T = waicNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:49:39,
## inlined from 'nimSmartPtr<T>& nimSmartPtr<T>::operator=(const nimSmartPtr<T>&) [with T = waicNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:79:11,
## inlined from 'virtual nimSmartPtr<waicNimbleList> waicClass::calculateWAIC(int, double)' at P_1_MCMC.cpp:1360:17:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '*<return-value>.nimSmartPtr<waicNimbleList>::realPtr' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'nimSmartPtr<T>::~nimSmartPtr() [with T = waicNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:116:29,
## inlined from 'virtual nimSmartPtr<waicNimbleList> waicClass::calculateWAIC(int, double)' at P_1_MCMC.cpp:1360:17:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '<unknown>' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## P_1_MCMC.cpp: In member function 'virtual nimSmartPtr<waicNimbleList> waicClass::calculateWAIC(int, double)':
## P_1_MCMC.cpp:1360:17: note: returned from 'void* operator new(std::size_t)'
## 1360 | Interm_84 = new waicNimbleList;
## | ^~~~~~~~~~~~~~
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'nimSmartPtr<T>::~nimSmartPtr() [with T = waicNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:116:29,
## inlined from 'SEXPREC* CALL_waicClass_get(SEXP)' at P_1_MCMC.cpp:1476:79:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '<anonymous>.nimSmartPtr<waicNimbleList>::realPtr' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'nimSmartPtr<T>::~nimSmartPtr() [with T = waicDetailsNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:116:29,
## inlined from 'SEXPREC* CALL_waicClass_getDetails(SEXP, SEXP)' at P_1_MCMC.cpp:1494:86:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '<anonymous>.nimSmartPtr<waicDetailsNimbleList>::realPtr' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'nimSmartPtr<T>::~nimSmartPtr() [with T = waicNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:116:29,
## inlined from 'SEXPREC* CALL_waicClass_calculateWAIC(SEXP, SEXP, SEXP)' at P_1_MCMC.cpp:1523:89:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '<anonymous>.nimSmartPtr<waicNimbleList>::realPtr' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'nimSmartPtr<T>::~nimSmartPtr() [with T = waicNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:116:29,
## inlined from 'double MCMC::calculateWAIC(int)' at P_1_MCMC.cpp:1724:37:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '<anonymous>.nimSmartPtr<waicNimbleList>::realPtr' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'nimSmartPtr<T>::~nimSmartPtr() [with T = waicNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:116:29,
## inlined from 'double MCMC::calculateWAIC(int)' at P_1_MCMC.cpp:1726:1:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '<anonymous>.nimSmartPtr<waicNimbleList>::realPtr' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'nimSmartPtr<T>::~nimSmartPtr() [with T = waicNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:116:29,
## inlined from 'nimSmartPtr<waicNimbleList> MCMC::getWAIC()' at P_1_MCMC.cpp:1731:31:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '<anonymous>.nimSmartPtr<waicNimbleList>::realPtr' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'nimSmartPtr<T>::~nimSmartPtr() [with T = waicNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:116:29,
## inlined from 'nimSmartPtr<waicNimbleList> MCMC::getWAIC()' at P_1_MCMC.cpp:1735:18:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '<unknown>' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## P_1_MCMC.cpp: In member function 'nimSmartPtr<waicNimbleList> MCMC::getWAIC()':
## P_1_MCMC.cpp:1735:18: note: returned from 'void* operator new(std::size_t)'
## 1735 | Interm_57 = new waicNimbleList;
## | ^~~~~~~~~~~~~~
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'nimSmartPtr<T>::~nimSmartPtr() [with T = waicNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:116:29,
## inlined from 'nimSmartPtr<waicNimbleList> MCMC::getWAIC()' at P_1_MCMC.cpp:1741:1:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '<unknown>' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## P_1_MCMC.cpp: In member function 'nimSmartPtr<waicNimbleList> MCMC::getWAIC()':
## P_1_MCMC.cpp:1735:18: note: returned from 'void* operator new(std::size_t)'
## 1735 | Interm_57 = new waicNimbleList;
## | ^~~~~~~~~~~~~~
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'nimSmartPtr<T>::~nimSmartPtr() [with T = waicNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:116:29,
## inlined from 'nimSmartPtr<waicNimbleList> MCMC::getWAIC()' at P_1_MCMC.cpp:1741:1:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '<anonymous>.nimSmartPtr<waicNimbleList>::realPtr' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'nimSmartPtr<T>::~nimSmartPtr() [with T = waicDetailsNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:116:29,
## inlined from 'nimSmartPtr<waicDetailsNimbleList> MCMC::getWAICdetails(bool)' at P_1_MCMC.cpp:1746:38:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '<anonymous>.nimSmartPtr<waicDetailsNimbleList>::realPtr' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'nimSmartPtr<T>::~nimSmartPtr() [with T = waicDetailsNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:116:29,
## inlined from 'nimSmartPtr<waicDetailsNimbleList> MCMC::getWAICdetails(bool)' at P_1_MCMC.cpp:1750:18:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '<unknown>' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## P_1_MCMC.cpp: In member function 'nimSmartPtr<waicDetailsNimbleList> MCMC::getWAICdetails(bool)':
## P_1_MCMC.cpp:1750:18: note: returned from 'void* operator new(std::size_t)'
## 1750 | Interm_59 = new waicDetailsNimbleList;
## | ^~~~~~~~~~~~~~~~~~~~~
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'nimSmartPtr<T>::~nimSmartPtr() [with T = waicDetailsNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:116:29,
## inlined from 'nimSmartPtr<waicDetailsNimbleList> MCMC::getWAICdetails(bool)' at P_1_MCMC.cpp:1758:1:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '<unknown>' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## P_1_MCMC.cpp: In member function 'nimSmartPtr<waicDetailsNimbleList> MCMC::getWAICdetails(bool)':
## P_1_MCMC.cpp:1750:18: note: returned from 'void* operator new(std::size_t)'
## 1750 | Interm_59 = new waicDetailsNimbleList;
## | ^~~~~~~~~~~~~~~~~~~~~
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'nimSmartPtr<T>::~nimSmartPtr() [with T = waicDetailsNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:116:29,
## inlined from 'nimSmartPtr<waicDetailsNimbleList> MCMC::getWAICdetails(bool)' at P_1_MCMC.cpp:1758:1:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '<anonymous>.nimSmartPtr<waicDetailsNimbleList>::realPtr' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'nimSmartPtr<T>::~nimSmartPtr() [with T = waicNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:116:29,
## inlined from 'SEXPREC* CALL_MCMC_getWAIC(SEXP)' at P_1_MCMC.cpp:1913:78:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '<anonymous>.nimSmartPtr<waicNimbleList>::realPtr' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'nimSmartPtr<T>::~nimSmartPtr() [with T = waicDetailsNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:116:29,
## inlined from 'SEXPREC* CALL_MCMC_getWAICdetails(SEXP, SEXP)' at P_1_MCMC.cpp:1931:85:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '<anonymous>.nimSmartPtr<waicDetailsNimbleList>::realPtr' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'void nimSmartPtr<T>::setPtr(const nimSmartPtr<T>&) [with T = waicNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:49:39,
## inlined from 'nimSmartPtr<T>& nimSmartPtr<T>::operator=(const nimSmartPtr<T>&) [with T = waicNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:79:11,
## inlined from 'virtual nimSmartPtr<waicNimbleList> waicClass::get()' at P_1_MCMC.cpp:1271:14:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '*<return-value>.nimSmartPtr<waicNimbleList>::realPtr' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'nimSmartPtr<T>::~nimSmartPtr() [with T = waicNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:116:29,
## inlined from 'virtual nimSmartPtr<waicNimbleList> waicClass::get()' at P_1_MCMC.cpp:1271:14:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '<unknown>' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## P_1_MCMC.cpp: In member function 'virtual nimSmartPtr<waicNimbleList> waicClass::get()':
## P_1_MCMC.cpp:1271:14: note: returned from 'void* operator new(std::size_t)'
## 1271 | output = new waicNimbleList;
## | ^~~~~~~~~~~~~~
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'void nimSmartPtr<T>::setPtr(const nimSmartPtr<T>&) [with T = waicDetailsNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:49:39,
## inlined from 'nimSmartPtr<T>& nimSmartPtr<T>::operator=(const nimSmartPtr<T>&) [with T = waicDetailsNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:79:11,
## inlined from 'virtual nimSmartPtr<waicDetailsNimbleList> waicClass::getDetails(bool)' at P_1_MCMC.cpp:1297:14:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '*<return-value>.nimSmartPtr<waicDetailsNimbleList>::realPtr' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## In destructor 'virtual pointedToBase::~pointedToBase()',
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:144:14,
## inlined from 'void pointedToBase::removeWatcher()' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:132:8,
## inlined from 'nimSmartPtr<T>::~nimSmartPtr() [with T = waicDetailsNimbleList]' at C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:116:29,
## inlined from 'virtual nimSmartPtr<waicDetailsNimbleList> waicClass::getDetails(bool)' at P_1_MCMC.cpp:1297:14:
## C:\Software\R\R-4.3.3\library\nimble\include/nimble/smartPtrs.h:147:29: warning: 'void operator delete(void*)' called on pointer '<unknown>' with nonzero offset 56 [-Wfree-nonheap-object]
## 147 | virtual ~pointedToBase() {};
## | ^
## P_1_MCMC.cpp: In member function 'virtual nimSmartPtr<waicDetailsNimbleList> waicClass::getDetails(bool)':
## P_1_MCMC.cpp:1297:14: note: returned from 'void* operator new(std::size_t)'
## 1297 | output = new waicDetailsNimbleList;
## | ^~~~~~~~~~~~~~~~~~~~~## running chain 1...
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|## running chain 2...
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## [Warning] There are 20 individual pWAIC values that are greater than 0.4. This may indicate that the WAIC estimate is unstable (Vehtari et al., 2017), at least in cases without grouping of data nodes or multivariate data nodes.nimbleOptions(showCompilerOutput = FALSE)
# See https://stats.stackexchange.com/q/304958 for discussion about WAIC warning
## View raw list output returned my the MCMC sampler ---------------------------
# Sampled values from the posterior distribution (n = n_iter - n_burin = 5,000)
utils::head(as.data.frame.matrix(bym_samples$samples$chain1), n = 10) # Summarised sample values from the posterior distribution (mean, median, std and quantiles 2.5 and 97.5)
utils::head(bym_samples$summary$all.chains, n = 100)## Mean Median St.Dev. 95%CI_low 95%CI_upp
## alpha 0.08888225 0.08930684 0.054984245 -0.02050222 0.1950861
## overall.sir 1.09460332 1.09341611 0.060121745 0.97970652 1.2154157
## residual.sir[1] 4.19434953 4.01978947 1.345700071 2.09085853 7.3304194
## residual.sir[2] 3.98422997 3.94304653 0.654064703 2.81293104 5.3804117
## residual.sir[3] 3.21591734 3.12435009 0.926864797 1.68128101 5.2882493
## residual.sir[4] 2.20131752 2.09530699 0.727440636 1.08713965 3.8855379
## residual.sir[5] 3.09086336 3.02883792 0.721573206 1.85486951 4.6868241
## residual.sir[6] 3.05292076 2.92076194 1.081827752 1.32380660 5.5151877
## residual.sir[7] 2.79297151 2.75102063 0.514265726 1.90344820 3.9129221
## residual.sir[8] 2.77959150 2.64308822 1.046786258 1.12588491 5.2024974
## residual.sir[9] 2.15076234 2.05961848 0.658081077 1.12995935 3.6847060
## residual.sir[10] 2.63240992 2.58951711 0.539666035 1.70527042 3.7995918
## residual.sir[11] 2.69481194 2.61806477 0.751047418 1.43471827 4.3582448
## residual.sir[12] 2.77096704 2.64319965 0.963928987 1.26321280 5.0158753
## residual.sir[13] 2.34322035 2.20479030 0.898523370 1.01142511 4.4612676
## residual.sir[14] 1.52416634 1.45119930 0.512304980 0.73588283 2.7167010
## residual.sir[15] 1.80033048 1.76178865 0.432697630 1.06819010 2.7594371
## residual.sir[16] 1.76537128 1.71575566 0.438046810 1.04536880 2.7533227
## residual.sir[17] 1.88321201 1.80073372 0.606796901 0.94866697 3.3128443
## residual.sir[18] 1.18518346 1.14546659 0.335929601 0.64852732 1.9523539
## residual.sir[19] 1.85583331 1.82003498 0.445506568 1.09196387 2.8321698
## residual.sir[20] 1.34513503 1.29242071 0.419249990 0.67763496 2.3088171
## residual.sir[21] 1.28260716 1.25641890 0.299746288 0.77646827 1.9408987
## residual.sir[22] 1.31154786 1.29716693 0.228629297 0.90402982 1.7988778
## residual.sir[23] 1.11280477 1.08553367 0.280295857 0.64678173 1.7345997
## residual.sir[24] 0.76302355 0.73788680 0.203904875 0.43731386 1.2314145
## residual.sir[25] 1.14798957 1.12903926 0.237490983 0.73873557 1.6702256
## residual.sir[26] 1.05891204 1.03586443 0.247566029 0.64165582 1.6093600
## residual.sir[27] 0.93797334 0.90417766 0.279611066 0.49251495 1.5850276
## residual.sir[28] 0.99618800 0.96980534 0.258330736 0.56853785 1.5732725
## residual.sir[29] 1.06857011 1.05301637 0.198565674 0.72443589 1.5012238
## residual.sir[30] 0.74626759 0.72547545 0.188672452 0.44023846 1.1762212
## residual.sir[31] 0.83967779 0.81038068 0.235647813 0.46328729 1.3811640
## residual.sir[32] 1.01496591 0.95714181 0.379752297 0.44818625 1.9124147
## residual.sir[33] 0.82437229 0.79598959 0.239200301 0.44011285 1.3689848
## residual.sir[34] 0.70683634 0.68989463 0.171421856 0.42291093 1.0937958
## residual.sir[35] 0.81201482 0.79082735 0.203357353 0.47172132 1.2606303
## residual.sir[36] 0.75176577 0.72723699 0.224712819 0.38582863 1.2644956
## residual.sir[37] 0.72799869 0.70845597 0.186846988 0.41760312 1.1492603
## residual.sir[38] 0.59830951 0.57774309 0.163896525 0.33726396 0.9804305
## residual.sir[39] 0.76688129 0.73701892 0.239843769 0.38282293 1.3130773
## residual.sir[40] 0.60095109 0.57228065 0.203947985 0.28711665 1.0798193
## residual.sir[41] 0.49746881 0.48621237 0.121687332 0.29293678 0.7689146
## residual.sir[42] 0.52911524 0.51670271 0.134269622 0.30345724 0.8257039
## residual.sir[43] 0.67956629 0.64804483 0.233171717 0.31647341 1.2208839
## residual.sir[44] 0.45326410 0.44085065 0.121084933 0.25350563 0.7255143
## residual.sir[45] 0.40860069 0.40315251 0.079934476 0.26717590 0.5792874
## residual.sir[46] 0.49794373 0.48161625 0.144892793 0.26163563 0.8270021
## residual.sir[47] 0.46101428 0.44275794 0.144409459 0.23308434 0.7917519
## residual.sir[48] 0.40765878 0.39054498 0.132278576 0.19737356 0.7127791
## residual.sir[49] 0.32920000 0.32607601 0.051890556 0.23629258 0.4398097
## residual.sir[50] 0.42971437 0.41710674 0.126353077 0.22114717 0.7136726
## residual.sir[51] 0.48466454 0.45741486 0.183549647 0.20904924 0.9168867
## residual.sir[52] 0.44294682 0.41708205 0.170492503 0.18867031 0.8454339
## residual.sir[53] 0.37186436 0.35312577 0.136078109 0.16213393 0.6852044
## residual.sir[54] 0.38469759 0.36893600 0.128440035 0.18019660 0.6785512
## residual.sir[55] 0.55345997 0.53296943 0.176834449 0.26779562 0.9560861
## residual.sir[56] 0.74947383 0.70815494 0.273826363 0.33876662 1.4008366
## residual.sir.gt1[1] 0.99995000 1.00000000 0.007070926 1.00000000 1.0000000
## residual.sir.gt1[2] 1.00000000 1.00000000 0.000000000 1.00000000 1.0000000
## residual.sir.gt1[3] 0.99951000 1.00000000 0.022130630 1.00000000 1.0000000
## residual.sir.gt1[4] 0.98601000 1.00000000 0.117449640 1.00000000 1.0000000
## residual.sir.gt1[5] 0.99998000 1.00000000 0.004472114 1.00000000 1.0000000
## residual.sir.gt1[6] 0.99510000 1.00000000 0.069828639 1.00000000 1.0000000
## residual.sir.gt1[7] 1.00000000 1.00000000 0.000000000 1.00000000 1.0000000
## residual.sir.gt1[8] 0.98578000 1.00000000 0.118397347 1.00000000 1.0000000
## residual.sir.gt1[9] 0.99032000 1.00000000 0.097910130 1.00000000 1.0000000
## residual.sir.gt1[10] 0.99999000 1.00000000 0.003162278 1.00000000 1.0000000
## residual.sir.gt1[11] 0.99848000 1.00000000 0.038957731 1.00000000 1.0000000
## residual.sir.gt1[12] 0.99354000 1.00000000 0.080114497 1.00000000 1.0000000
## residual.sir.gt1[13] 0.97657000 1.00000000 0.151265541 1.00000000 1.0000000
## residual.sir.gt1[14] 0.86383000 1.00000000 0.342970709 0.00000000 1.0000000
## residual.sir.gt1[15] 0.98676000 1.00000000 0.114301501 1.00000000 1.0000000
## residual.sir.gt1[16] 0.98287000 1.00000000 0.129756431 1.00000000 1.0000000
## residual.sir.gt1[17] 0.96501000 1.00000000 0.183755374 0.00000000 1.0000000
## residual.sir.gt1[18] 0.68291000 1.00000000 0.465345138 0.00000000 1.0000000
## residual.sir.gt1[19] 0.98864000 1.00000000 0.105976708 1.00000000 1.0000000
## residual.sir.gt1[20] 0.79059000 1.00000000 0.406889552 0.00000000 1.0000000
## residual.sir.gt1[21] 0.82857000 1.00000000 0.376886157 0.00000000 1.0000000
## residual.sir.gt1[22] 0.92441000 1.00000000 0.264342298 0.00000000 1.0000000
## residual.sir.gt1[23] 0.62747000 1.00000000 0.483480854 0.00000000 1.0000000
## residual.sir.gt1[24] 0.12165000 0.00000000 0.326882771 0.00000000 1.0000000
## residual.sir.gt1[25] 0.71812000 1.00000000 0.449917426 0.00000000 1.0000000
## residual.sir.gt1[26] 0.55982000 1.00000000 0.496411152 0.00000000 1.0000000
## residual.sir.gt1[27] 0.36594000 0.00000000 0.481695170 0.00000000 1.0000000
## residual.sir.gt1[28] 0.45324000 0.00000000 0.497811190 0.00000000 1.0000000
## residual.sir.gt1[29] 0.60805000 1.00000000 0.488188059 0.00000000 1.0000000
## residual.sir.gt1[30] 0.09557000 0.00000000 0.294002108 0.00000000 1.0000000
## residual.sir.gt1[31] 0.22206000 0.00000000 0.415633353 0.00000000 1.0000000
## residual.sir.gt1[32] 0.45396000 0.00000000 0.497878296 0.00000000 1.0000000
## residual.sir.gt1[33] 0.21048000 0.00000000 0.407651605 0.00000000 1.0000000
## residual.sir.gt1[34] 0.05750000 0.00000000 0.232796675 0.00000000 1.0000000
## residual.sir.gt1[35] 0.17176000 0.00000000 0.377173601 0.00000000 1.0000000
## residual.sir.gt1[36] 0.13281000 0.00000000 0.339370971 0.00000000 1.0000000
## residual.sir.gt1[37] 0.08322000 0.00000000 0.276215848 0.00000000 1.0000000
## residual.sir.gt1[38] 0.02094000 0.00000000 0.143184222 0.00000000 0.0000000
## residual.sir.gt1[39] 0.15846000 0.00000000 0.365173605 0.00000000 1.0000000
## residual.sir.gt1[40] 0.04188000 0.00000000 0.200315918 0.00000000 1.0000000
## residual.sir.gt1[41] 0.00070000 0.00000000 0.026448384 0.00000000 0.0000000
## residual.sir.gt1[42] 0.00259000 0.00000000 0.050826349 0.00000000 0.0000000### Check model convergence using the Gelman-Rubin convergence diagnostic ------
# The Gelman-Rubin compares the within chain variation to the between chain
# variation for each parameter. If the chains have converged on the stationary
# distribution, then it should be impossible to distinguish between the chains.
# It is a ratio measures, and values < 1.1 are considered evidence of chain
# convergence.
# gr.diag <- coda::gelman.diag(bym_samples$samples, multivariate = FALSE)
#
# # List all parameters that are not posterior probabilities
#
# gr.diag.inc <- stringr::str_subset(base::dimnames(gr.diag$psrf)[[1]], "gt1|lt1", negate = TRUE)
#
# ## Are all values < 1.1? If yes, this is evidence of convergence ---------------
#
# base::all(gr.diag$psrf[gr.diag.inc,"Point est."] < 1.1) # Don't check posterior probabilities - these are threshold nodes and do not need to converge
#
# ## If no, then Which parameters have a value >= 1.1? ---------------------------
#
# base::which(gr.diag$psrf[gr.diag.inc,"Point est."] > 1.1)
### Check trace, density and autocorrelation plots for convergence -------------
# Import MCMC samples into a ggs that can be used with ggs_* graphical functions
bym_ggmcmc <- ggmcmc::ggs(bym_samples$samples)
## Check the overall SIR -------------------------------------------------------
# Generate trace, density and autocorrelation plots
trace_overall.sir <- bym_ggmcmc %>%
dplyr::filter(Parameter == "overall.sir") %>%
ggmcmc::ggs_traceplot() +
ggplot2::theme_bw() +
ggplot2::theme(text = ggplot2::element_text(size = 12))
density_overall.sir <- bym_ggmcmc %>%
dplyr::filter(Parameter == "overall.sir") %>%
ggmcmc::ggs_density() +
ggplot2::theme_bw() +
ggplot2::theme(text = ggplot2::element_text(size = 12))
acf_overall.sir <- bym_ggmcmc %>%
dplyr::filter(Parameter == "overall.sir") %>%
ggmcmc::ggs_autocorrelation() +
ggplot2::theme_bw() +
ggplot2::theme(text = element_text(size = 12))
# Plot visual convergence diagnostics
base::print((trace_overall.sir | density_overall.sir) / acf_overall.sir) +
patchwork::plot_annotation(title = "Trace, density and autocorrelation plots for overall SIR",
theme = ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)))## Check the spatially variance component --------------------------------------
# Generate trace, density and autocorrelation plots
trace_sigma2.s <- bym_ggmcmc %>%
dplyr::filter(Parameter == "sigma2.s") %>%
ggmcmc::ggs_traceplot() +
ggplot2::theme_bw() +
ggplot2::theme(text = ggplot2::element_text(size = 12))
density_sigma2.s <- bym_ggmcmc %>%
dplyr::filter(Parameter == "sigma2.s") %>%
ggmcmc::ggs_density() +
ggplot2::theme_bw() +
ggplot2::theme(text = ggplot2::element_text(size = 12))
acf_sigma2.s <- bym_ggmcmc %>%
filter(Parameter == "sigma2.s") %>%
ggmcmc::ggs_autocorrelation() +
ggplot2::theme_bw() +
ggplot2::theme(text = ggplot2::element_text(size = 12))
# Plot visual convergence diagnostics
base::print((trace_sigma2.s | density_sigma2.s)/acf_sigma2.s) +
patchwork::plot_annotation(title = "Trace, density and autocorrelation plots for spatially correlated variance component (s)",
theme = ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)))## Check the uncorrelated variance component -----------------------------------
# Generate trace, density and autocorrelation plots
trace_sigma2.u <- bym_ggmcmc %>%
dplyr::filter(Parameter == "sigma2.u") %>%
ggmcmc::ggs_traceplot() +
ggplot2::theme_bw() +
ggplot2::theme(text = ggplot2::element_text(size = 12))
density_sigma2.u <- bym_ggmcmc %>%
dplyr::filter(Parameter == "sigma2.u") %>%
ggmcmc::ggs_density() +
ggplot2::theme_bw() +
ggplot2::theme(text = ggplot2::element_text(size = 12))
acf_sigma2.u <- bym_ggmcmc %>%
filter(Parameter == "sigma2.u") %>%
ggmcmc::ggs_autocorrelation() +
ggplot2::theme_bw() +
ggplot2::theme(text = ggplot2::element_text(size = 12))
base::print((trace_sigma2.u | density_sigma2.u) /acf_sigma2.u) +
patchwork::plot_annotation(title = "Trace, density and autocorrelation plots for uncorrelated variance component (u)",
theme = ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)))## Check the smoothed SIR for Perth-Kinross (29) -------------------------------
# This district has the most neighbours
# Generate trace, density and autocorrelation plots
trace_smoothed.sir_29 <- bym_ggmcmc %>%
dplyr::filter(Parameter == "smoothed.sir[29]") %>%
ggmcmc::ggs_traceplot() +
ggplot2::theme_bw() +
ggplot2::theme(text = ggplot2::element_text(size = 12))
density_smoothed.sir_29 <- bym_ggmcmc %>%
dplyr::filter(Parameter == "smoothed.sir[29]") %>%
ggmcmc::ggs_density() +
ggplot2::theme_bw() +
ggplot2::theme(text = ggplot2::element_text(size = 12))
acf_smoothed.sir_29 <- bym_ggmcmc %>%
filter(Parameter == "smoothed.sir[29]") %>%
ggmcmc::ggs_autocorrelation() +
ggplot2::theme_bw() +
ggplot2::theme(text = ggplot2::element_text(size = 12))
# Plot visual convergence diagnostics
base::print((trace_smoothed.sir_29 | density_smoothed.sir_29) / acf_smoothed.sir_29) +
patchwork::plot_annotation(title = "Trace, density and autocorrelation plots for Perth-Kinross smoothed SIR",
theme = ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)))### Summarise across MCMC samples to obtained BYM model estimates --------------
## Extract alpha, sigma2.s and sigma2.u ----------------------------------------
bym_ests <- bym_samples$summary$all.chains[c("alpha",
"sigma2.s",
"sigma2.u"),
c("Median",
"95%CI_low",
"95%CI_upp")] %>%
# Convert matrix to data frame
base::as.data.frame() %>%
# Rename columns
stats::setNames(c("median", "lower_95_ci", "upper_95_ci")) %>%
# Make the data frame row names an explicit variable called *parameter*
tibble::rownames_to_column(var = "parameter") %>%
# Exponentiate the parameter and 95% CI estimates for alpha to obtain overall SIR effect measure
dplyr::mutate(dplyr::across(c(median, lower_95_ci, upper_95_ci), ~ base::ifelse(parameter == "alpha", base::exp(.), NA), .names = "{.col}_exp")) %>%
# Format estimates as text variables with three decimal places
dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ base::ifelse(base::is.na(.), "", base::format(base::round(., 3), nsmall = 3))))
### Print model estimates and WIAC ---------------------------------------------
base::print(list("Model estimate" = bym_ests,
"WAIC" = bym_samples$WAIC$WAIC))## $`Model estimate`
## parameter median lower_95_ci upper_95_ci median_exp lower_95_ci_exp upper_95_ci_exp
## 1 alpha 0.089 -0.021 0.195 1.093 0.980 1.215
## 2 sigma2.s 0.605 0.316 1.132
## 3 sigma2.u 0.011 0.003 0.066
##
## $WAIC
## [1] 293.8658### Summarise across MCMC samples to obtain smoothed SIR estimates -------------
# Create data set of smoothed and residual risk of lip cancer
bym_sf <- lip_sf %>%
# Smoothed SIR estimate
dplyr::bind_cols(bym_samples$summary$all.chains[base::paste0("smoothed.sir[", 1:N, "]"), c("Median", "95%CI_low", "95%CI_upp")] %>%
base::as.data.frame() %>%
stats::setNames(c("bym_smoothed_sir_est", "bym_smoothed_sir_l95", "bym_smoothed_sir_uci"))) %>%
# Residual SIR estimate
dplyr::bind_cols(bym_samples$summary$all.chains[base::paste0("residual.sir[", 1:N, "]"), c("Median", "95%CI_low", "95%CI_upp")] %>%
base::as.data.frame() %>%
stats::setNames(c("bym_residual_sir_est", "bym_residual_sir_l95", "bym_residual_sir_uci"))) %>%
# Residual s SIR estimate
dplyr::bind_cols(bym_samples$summary$all.chains[base::paste0("residual.sir.s[", 1:N, "]"), c("Median", "95%CI_low", "95%CI_upp")] %>%
base::as.data.frame() %>%
stats::setNames(c("bym_residual_sir_s_est", "bym_residual_sir_s_l95", "bym_residual_sir_s_uci"))) %>%
# Residual u SIR estimate
dplyr::bind_cols(bym_samples$summary$all.chains[base::paste0("residual.sir.u[", 1:N, "]"), c("Median", "95%CI_low", "95%CI_upp")] %>%
base::as.data.frame() %>%
stats::setNames(c("bym_residual_sir_u_est", "bym_residual_sir_u_l95", "bym_residual_sir_u_uci"))) %>%
# Posterior probabilities for residual risk > 1
dplyr::bind_cols(bym_samples$summary$all.chains[base::paste0("residual.sir.gt1[", 1:N, "]"), "Mean"] %>%
base::as.data.frame() %>%
stats::setNames("bym_residual_sir_gt1")) %>%
# Posterior probabilities for residual risk < 1
dplyr::bind_cols(bym_samples$summary$all.chains[base::paste0("residual.sir.lt1[", 1:N, "]"), "Mean"] %>%
base::as.data.frame() %>%
stats::setNames("bym_residual_sir_lt1")) %>%
# Annotate exceedence probabilities
dplyr::mutate(bym_residual_sir_cat = dplyr::case_when(bym_residual_sir_gt1 > 0.995 ~ "++",
bym_residual_sir_gt1 > 0.975 ~ "+",
bym_residual_sir_lt1 > 0.995 ~ "--",
bym_residual_sir_lt1 > 0.975 ~ "-",
.default = ""))
### Map BYM smoothed SIR for districts -----------------------------------------
map_bym_residual_sir <- ggplot() +
# Add SIR estimates
geom_sf(data = bym_sf,
aes(fill = bym_residual_sir_est),
col = "grey50",
linewidth = 0.1) +
# Add SIR posterior probabilities
ggplot2::geom_point(data = bym_sf %>% dplyr::filter(bym_residual_sir_cat != ""),
mapping = ggplot2::aes(shape = bym_residual_sir_cat,
geometry = geom),
fill = "black",
stat = "sf_coordinates") +
# Add Scotland boundary
ggplot2::geom_sf(data = sct_sf,
fill = NA,
colour = "Black",
linewidth = 1) +
# Add fill aesthetic gradient mapping
ggplot2::scale_fill_gradient2("Smoothed SIR",
low = col_pal["low"], # blue
mid = col_pal["mid"], # yellow
high = col_pal["high"], # darkblue
midpoint = 0,
na.value = NA,
limits = c(base::exp(-max_log_sir), base::exp(max_log_sir)),
breaks = base::exp(base::seq(-max_log_sir, max_log_sir, 0.5)),
labels = base::format(base::round(exp(seq(-max_log_sir, max_log_sir, 0.5)), 2), nsmall = 2),
trans = "log") +
# Add shape aesthic mapping
ggplot2::scale_shape_manual("Posterior probabilities",
values = c("++" = 24,
"+" = 2,
"-" = 6,
"--" = 25),
breaks = c("++", "+", "-", "--"),
labels = c("++" = "SIR > 1, p > 97.5%",
"+" = "SIR > 1, p > 95.0%",
"-" = "SIR < 1, p > 95.0%",
"--" = "SIR < 1, p > 97.5%")) +
# Format legends
guides(fill = guide_colourbar(frame.colour = "black",
frame.linewidth = 0.25,
ticks.colour = "black",
ticks.linewidth = 0.25,
order = 1),
shape = guide_legend(order = 2)) +
# Format map
theme_bw() +
theme(axis.title = ggplot2::element_blank(),
axis.text = ggplot2::element_blank(),
axis.ticks = ggplot2::element_blank(),
legend.position = "inside",
legend.margin = ggplot2::margin(0, 0, 0, 0),
legend.position.inside = c(0.15, 0.88),
panel.grid = ggplot2::element_blank(),
text = ggplot2::element_text(size = 12))
## Output observed and BYM smoothed SIR for lip cancer risk --------------------
(map_obs_sir | map_bym_residual_sir) +
patchwork::plot_annotation(title = "Observed and BYM smoothed SIR for male lip cancer in Scotland (1975-1980)",
theme = theme(plot.title = element_text(hjust = 0.5)))Clayton D, Kaldor J. Empirical Bayes estimates of age standardised relative risks for use in disease mapping. Biometrics 1987;43:671–681. doi: 10.2307/2532003
Cressie NAC. Statistics for spatial data. New York: John Wiley & Sons; 1993. p. 537 (Table 7.2). doi: 10.1002/9781119115151
Stern H, Cressie N. Inferences for extremes in disease mapping. In Lawson AB, Biggeri A, Böhning D, et al., editors. Disease mapping and risk assessment for public health. New York: Wiley; 1999. pp. 68–69 (Table 5.1). Available at Storage General (614.42 41)